Coral Performance Following Reciprocal Transplant

Setup

Load Packages

##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("cowplot")) install.packages("cowplot")
if (!require("lubridate")) install.packages("lubridate")
if (!require("corrplot")) install.packages("corrplot")
if (!require("lme4")) install.packages("lme4")
if (!require("vegan")) install.packages("vegan")
if (!require("DHARMa")) install.packages("DHARMa")
if (!require("effectsize")) install.packages("effectsize")
if (!require("emmeans")) install.packages("emmeans")
if (!require("lmerTest")) install.packages("lmerTest")
if (!require("tidyr")) install.packages("tidyr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("plotrix")) install.packages("plotrix")
if (!require("Rmisc")) install.packages("Rmisc")
if (!require("ggpubr")) install.packages("ggpubr")

##Load Packages
library(ggplot2) #Required for ggplots
library(cowplot) #Required for plotting panel figures
library(lubridate) #Required for date and time formatting
library(corrplot) #Required for correlation plot
library(lme4) #Required for mixed effects modeling
library(vegan) #Required for multivariate analysis PERM
library(DHARMa) #Required to check residuals of mixed effects models
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons 
library(lmerTest) #Required for p values with lme4 model summaries
library(tidyr) #Required for pivot_wider function
library(dplyr) #Required for group_by function
library(plotrix) #Required for Standard error function
library(Rmisc) #Required for SummarySE function
library(ggpubr) #Required for stat brackets

Graphing Parameters

#Note: Run "Graphing Parameters" section from 01_ExperimentalSetup.Rmd file

Sample Data and Metadata

Load and Organize Data

##Growth
Length<-read.csv("Data/LengthFull.csv", header=TRUE)
Dates<-read.csv("Data/BonaireDates.csv", header=TRUE)

##Set date variables
Dates$Date<-as.Date(Dates$Date, format='%m/%d/%Y')

##Set factor variables
Length$Site<-factor(Length$Site, levels=c("KL", "SS"))
Length$Genotype<-factor(Length$Genotype, levels=c("AC8", "AC10", "AC12"))
Length$Orig<-factor(Length$Orig, levels=c("N", "T"))
Length$Origin<-factor(Length$Origin, levels=c("Native", "Transplant"))

##Add a Sample Set Variable
Length$Set<-paste(Length$Site, Length$Genotype, Length$Orig, sep=".")
Length$Set<-factor(Length$Set)

##Add Site.Orig variable
Length$Site.Orig<-paste(Length$Site, Length$Orig, sep=".")
Length$Site.Orig<-factor(Length$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))


##Bleaching Metrics
#Note: Physiological metrics calculated in 02_PhysiologyMetrics.R file
Thermal<-read.csv("Outputs/ThermData.csv", header=TRUE)

##Set factor variables
Thermal$TimeP<-factor(Thermal$TimeP, levels=c("IN", "TP1", "TP2", "TP3", "TP4"))
Thermal$Site<-factor(Thermal$Site, levels=c("KL", "SS"))
Thermal$Genotype<-factor(Thermal$Genotype, levels=c("AC8", "AC10", "AC12"))
Thermal$Orig<-factor(Thermal$Orig, levels=c("N", "T"))
Thermal$Origin<-factor(Thermal$Origin, levels=c("Native", "Transplant"))

Growth

Growth Rate TP1 to TP2

Calculate Days

##Subset Dates for TP 1 and TP 2 Growth
Growth.Dates_TP1.2<-subset(Dates, Task=="Growth" & TimeP=="TP1" |
                       Task=="Growth" & TimeP=="TP2" )

##Convert to Wide format
Growth.Dates_TP1.2<-Growth.Dates_TP1.2 %>%  pivot_wider(names_from = TimeP, values_from = Date)

##Calculate Number of Days
#TP 1 to TP2
Growth.Dates_TP1.2$TP1.2_days<-time_length(Growth.Dates_TP1.2$TP2-Growth.Dates_TP1.2$TP1, unit="days")

Calculate Total Linear Extension

##Merge Length Data with Growth Dates
#Merges by Site column
#Adds TP1.2_days column
Growth_TP1.2<-merge(Length, Growth.Dates_TP1.2[,c(-2)], all.x=TRUE)

##Calculate Linear Extension (cm)
Growth_TP1.2$Ext_cm<-Growth_TP1.2$TL.TP2_cm-Growth_TP1.2$TL.TP1_cm

##Calculate Growth_TP1.2 Rate (cm/day)
Growth_TP1.2$TLE_cm.day<-Growth_TP1.2$Ext_cm/Growth_TP1.2$TP1.2_days

Data QC

##Remove Corals not Measured in T1 to T2
Growth_TP1.2<-Growth_TP1.2 %>% drop_na(TLE_cm.day)

##Sample Size per Set
Growth_TP1.2 %>% 
  dplyr::group_by(Set) %>%
  dplyr::summarize(SampleSize = length(Set))

##Check for Outliers
ggplot(Growth_TP1.2, aes(x=Set, y=TLE_cm.day)) + 
  geom_boxplot(alpha=0.5, shape=2, outlier.shape = NA)+
  geom_jitter(shape=16, position=position_jitter(0.1))+
  theme(axis.text.x = element_text(angle = 90))

Potential outliers have been re-measured and therefore will be retained. n=10 per Site, Genotype, Origin group

Compare Growth Rate

#Check normality
hist(Growth_TP1.2$TLE_cm.day)

shapiro.test(Growth_TP1.2$TLE_cm.day)

    Shapiro-Wilk normality test

data:  Growth_TP1.2$TLE_cm.day
W = 0.92803, p-value = 7.262e-06

#Not normal

hist(log(Growth_TP1.2$TLE_cm.day+1))

shapiro.test(log(Growth_TP1.2$TLE_cm.day+1))

    Shapiro-Wilk normality test

data:  log(Growth_TP1.2$TLE_cm.day + 1)
W = 0.94335, p-value = 7.248e-05

##Still not normal

##Compare generalized linear mixed effects models
TLE_TP1.2.glmr.gaus<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=gaussian(link="identity"))

Warning: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly

TLE_TP1.2.glmr.gaus.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=gaussian(link="log"))
TLE_TP1.2.glmr.gaus.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=gaussian(link="inverse"))

TLE_TP1.2.glmr.gam.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=Gamma(link="log"))
TLE_TP1.2.glmr.gam.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=Gamma(link="inverse"))

TLE_TP1.2.glmr.inv<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=inverse.gaussian(link="1/mu^2"))

boundary (singular) fit: see help('isSingular')

TLE_TP1.2.glmr.inv.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=inverse.gaussian(link="log"))
TLE_TP1.2.glmr.inv.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP1.2, family=inverse.gaussian(link="inverse"))

AIC(TLE_TP1.2.glmr.gaus, TLE_TP1.2.glmr.gaus.l, TLE_TP1.2.glmr.gaus.i, TLE_TP1.2.glmr.gam.l, TLE_TP1.2.glmr.gam.i, TLE_TP1.2.glmr.inv, TLE_TP1.2.glmr.inv.l, TLE_TP1.2.glmr.inv.i)

Gamma distribution with log-link has the lowest AIC.

##Check residuals
TLE_TP1.2.glmr.gam.l_res <- simulateResiduals(fittedModel = TLE_TP1.2.glmr.gam.l, plot = F)
plot(TLE_TP1.2.glmr.gam.l_res)

Compare residuals across models

##Check residuals
TLE_TP1.2.glmr.gaus.l_res <- simulateResiduals(fittedModel = TLE_TP1.2.glmr.gaus.l, plot = F)
plot(TLE_TP1.2.glmr.gaus.l_res)

Compare with log+1 transformed lmer model

LMER Model

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
TLE_TP1.2.lme<-lmer(log(TLE_cm.day+1)~Origin*Site+(1|Genotype), data=Growth_TP1.2)

##Check residuals
TLE_TP1.2.lme_res <- simulateResiduals(fittedModel = TLE_TP1.2.lme, plot = F)
plot(TLE_TP1.2.lme_res)

Better residuals with the log+1 transformed lmer model.

##Model results
summary(TLE_TP1.2.lme)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin * Site + (1 | Genotype)
   Data: Growth_TP1.2

REML criterion at convergence: -233.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0026 -0.5991 -0.1118  0.5959  2.4282 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002965 0.05445 
 Residual             0.006632 0.08144 
Number of obs: 120, groups:  Genotype, 3

Fixed effects:
                          Estimate Std. Error         df t value Pr(>|t|)   
(Intercept)               0.125179   0.034778   2.684719   3.599  0.04400 * 
OriginTransplant         -0.007538   0.021027 114.000000  -0.359  0.72063   
SiteSS                    0.062096   0.021027 114.000000   2.953  0.00382 **
OriginTransplant:SiteSS  -0.035167   0.029737 114.000000  -1.183  0.23943   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.302              
SiteSS      -0.302  0.500       
OrgnTrn:SSS  0.214 -0.707 -0.707

eta_squared(TLE_TP1.2.lme)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.02 | [0.00, 1.00]
Site        |           0.07 | [0.01, 1.00]
Origin:Site |           0.01 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
TLE_TP1.2.res<-data.frame(summary(TLE_TP1.2.lme)$coefficients[-1,])
TLE_TP1.2.res$Predictor<-c("Origin", "Site", "Origin x Site")
TLE_TP1.2.res$EtaSq<-c(eta_squared(TLE_TP1.2.lme)$Eta2)
TLE_TP1.2.res$Response<-rep("Growth", nrow(TLE_TP1.2.res))
TLE_TP1.2.res$Timepoint<-rep("TP1_2", nrow(TLE_TP1.2.res))

Pairwise

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
TLE.pair_TP1.2<-emmeans(TLE_TP1.2.lme, pairwise~Origin | Site)
TLE.pair_TP1.2

$emmeans
Site = KL:
 Origin     emmean     SE   df  lower.CL upper.CL
 Native      0.125 0.0348 2.68  0.006770    0.244
 Transplant  0.118 0.0348 2.68 -0.000768    0.236

Site = SS:
 Origin     emmean     SE   df  lower.CL upper.CL
 Native      0.187 0.0348 2.68  0.068866    0.306
 Transplant  0.145 0.0348 2.68  0.026161    0.263

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate    SE  df t.ratio p.value
 Native - Transplant  0.00754 0.021 114   0.359  0.7206

Site = SS:
 contrast            estimate    SE  df t.ratio p.value
 Native - Transplant  0.04271 0.021 114   2.031  0.0446

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
TLE.pair.es_TP1.2<-data.frame(eff_size(TLE.pair_TP1.2, sigma=sigma(TLE_TP1.2.lme), edf=df.residual(TLE_TP1.2.lme)))

Since 'object' is a list, we are using the contrasts already present.

TLE.pair.es_TP1.2

TLE.pair.es_TP1.2<-TLE.pair.es_TP1.2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
TLE.pair_TP1.2.res<-merge(data.frame(TLE.pair_TP1.2$contrasts), TLE.pair.es_TP1.2[,-c(1)])
TLE.pair_TP1.2.res$Response<-rep("Growth", nrow(TLE.pair_TP1.2.res))
TLE.pair_TP1.2.res$TimeP<-rep("TP1_2", nrow(TLE.pair_TP1.2.res))

##Add Significance levels
TLE.pair_TP1.2.res$Sig<-ifelse(TLE.pair_TP1.2.res$p.value<0.001, "***",
                              ifelse(TLE.pair_TP1.2.res$p.value<0.01, "**",
                                     ifelse(TLE.pair_TP1.2.res$p.value<0.05, "*",
                                            ifelse(TLE.pair_TP1.2.res$p.value<0.1, "-", NA))))
TLE.pair_TP1.2.res

Plot Growth TP 1 to 2

##Summary statistics by Site and Origin
TP1.2_Growth.sum<-summarySE(Growth_TP1.2, measurevar="TLE_cm.day", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Growth Rate across Treatments
TP1.2_Growth.plot<-ggplot(TP1.2_Growth.sum, aes(x=Site, y=TLE_cm.day, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=TLE_cm.day-se, ymax=TLE_cm.day+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y=expression(paste('Growth Rate (cm day'^-1*")")), colour="Origin")+
  ylim(0, 0.25)+
  annotate("text", x=2, y=0.24, label="*", size=sig.sz, fontface="bold")+
  geom_bracket(xmin=1, xmax=2, y.position=0.01, label.size=levels.sz, label="**", inherit.aes = FALSE); TP1.2_Growth.plot

Percent Difference in Growth Rates

##Summary statistics by Genotype, Origin, and Site
TP1.2_Growth.sum<-summarySE(Growth_TP1.2, measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin", "Site"), na.rm=TRUE)

##Calculate Percent Difference between Native and Transplant by Genotype
TP1.2_Growth.dif<-TP1.2_Growth.sum[which(TP1.2_Growth.sum$Origin=="Native"), c(1,3,5)]
TP1.2_Growth.dif<-TP1.2_Growth.dif %>% dplyr::rename(TLE_N = TLE_cm.day)
TP1.2_Growth.dif$TLE_T<-TP1.2_Growth.sum$TLE_cm.day[which(TP1.2_Growth.sum$Origin=="Transplant")]

##SS
TP1.2_Growth.dif_SS<-subset(TP1.2_Growth.dif, Site=="SS")
TP1.2_Growth.dif_SS$Perc.Inc<-((TP1.2_Growth.dif_SS$TLE_N-TP1.2_Growth.dif_SS$TLE_T)/TP1.2_Growth.dif_SS$TLE_T)*100
TP1.2_Growth.dif_SS

##Mean and Standard Error
mean(TP1.2_Growth.dif_SS$Perc.Inc)

[1] 47.19841

std.error(TP1.2_Growth.dif_SS$Perc.Inc)

[1] 25.89777

##KL
TP1.2_Growth.dif_KL<-subset(TP1.2_Growth.dif, Site=="KL")
TP1.2_Growth.dif_KL$Perc.Inc<-((TP1.2_Growth.dif_KL$TLE_N-TP1.2_Growth.dif_KL$TLE_T)/TP1.2_Growth.dif_KL$TLE_T)*100
TP1.2_Growth.dif_KL

##Mean and Standard Error
mean(TP1.2_Growth.dif_KL$Perc.Inc)

[1] 4.187086

std.error(TP1.2_Growth.dif_KL$Perc.Inc)

[1] 8.336196

Growth Rate TP3 to TP4

Calculate Days

##Subset Dates for TP 3 and TP 4 Growth
Growth.Dates_TP3.4<-subset(Dates, Task=="Growth" & TimeP=="TP3" |
                       Task=="Growth" & TimeP=="TP4" )

##Convert to Wide format
Growth.Dates_TP3.4<-Growth.Dates_TP3.4 %>%  pivot_wider(names_from = TimeP, values_from = Date)

##Calculate Number of Days
#TP 3 to TP4
Growth.Dates_TP3.4$TP3.4_days<-time_length(Growth.Dates_TP3.4$TP4-Growth.Dates_TP3.4$TP3, unit="days")

Calculate Total Linear Extension

##Merge Length Data with Growth Dates
#Merges by Site column
#Adds TP3.4_days column
Growth_TP3.4<-merge(Length, Growth.Dates_TP3.4[,c(-2)], all.x=TRUE)

##Calculate Linear Extension (cm)
Growth_TP3.4$Ext_cm<-Growth_TP3.4$TL.TP4_cm-Growth_TP3.4$TL.TP3_cm

##Calculate Growth_TP3.4 Rate (cm/day)
Growth_TP3.4$TLE_cm.day<-Growth_TP3.4$Ext_cm/Growth_TP3.4$TP3.4_days

Data QC

##Remove Corals not Measured in T1 to T2
Growth_TP3.4<-Growth_TP3.4 %>% drop_na(TLE_cm.day)

##Sample Size per Set
Growth_TP3.4 %>% 
  dplyr::group_by(Set) %>%
  dplyr::summarize(SampleSize = length(Set))

##Check for Outliers
ggplot(Growth_TP3.4, aes(x=Set, y=TLE_cm.day)) + 
  geom_boxplot(alpha=0.5, shape=2, outlier.shape = NA)+
  geom_jitter(shape=16, position=position_jitter(0.1))+
  theme(axis.text.x = element_text(angle = 90))

Compare Growth Rate

#Check normality
hist(Growth_TP3.4$TLE_cm.day)

shapiro.test(Growth_TP3.4$TLE_cm.day)

    Shapiro-Wilk normality test

data:  Growth_TP3.4$TLE_cm.day
W = 0.92767, p-value = 6.904e-06

#Not normal

hist(log(Growth_TP3.4$TLE_cm.day+1))

shapiro.test(log(Growth_TP3.4$TLE_cm.day+1))

    Shapiro-Wilk normality test

data:  log(Growth_TP3.4$TLE_cm.day + 1)
W = 0.94719, p-value = 0.0001351

##Still not normal

##Compare generalized linear mixed effects models
TLE_TP3.4.glmr.gaus<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=gaussian(link="identity"))

Warning: calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly

TLE_TP3.4.glmr.gaus.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=gaussian(link="log"))
TLE_TP3.4.glmr.gaus.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=gaussian(link="inverse"))

TLE_TP3.4.glmr.gam.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=Gamma(link="log"))
TLE_TP3.4.glmr.gam.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=Gamma(link="inverse"))

TLE_TP3.4.glmr.inv<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=inverse.gaussian(link="1/mu^2"))
TLE_TP3.4.glmr.inv.l<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=inverse.gaussian(link="log"))
TLE_TP3.4.glmr.inv.i<-glmer(TLE_cm.day~Origin*Site+(1|Genotype), data=Growth_TP3.4, family=inverse.gaussian(link="inverse"))

AIC(TLE_TP3.4.glmr.gaus, TLE_TP3.4.glmr.gaus.l, TLE_TP3.4.glmr.gaus.i, TLE_TP3.4.glmr.gam.l, TLE_TP3.4.glmr.gam.i, TLE_TP3.4.glmr.inv, TLE_TP3.4.glmr.inv.l, TLE_TP3.4.glmr.inv.i)

Gamma distribution with log-link has the lowest AIC.

##Check residuals
TLE_TP3.4.glmr.gam.l_res <- simulateResiduals(fittedModel = TLE_TP3.4.glmr.gam.l, plot = F)
plot(TLE_TP3.4.glmr.gam.l_res)

Compare residuals across models

##Check residuals
TLE_TP3.4.glmr.gaus.l_res <- simulateResiduals(fittedModel = TLE_TP3.4.glmr.gaus.l, plot = F)
plot(TLE_TP3.4.glmr.gaus.l_res)

Compare with log+1 transformed lmer model

LMER Model

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
TLE_TP3.4.lme<-lmer(log(TLE_cm.day+1)~Origin*Site+(1|Genotype), data=Growth_TP3.4)

##Check residuals
TLE_TP3.4.lme_res <- simulateResiduals(fittedModel = TLE_TP3.4.lme, plot = F)
plot(TLE_TP3.4.lme_res)

Better residuals with the log+1 transformed lmer model.

##Model results
summary(TLE_TP3.4.lme)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(TLE_cm.day + 1) ~ Origin * Site + (1 | Genotype)
   Data: Growth_TP3.4

REML criterion at convergence: -341.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.27432 -0.60263 -0.01933  0.56722  2.45603 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002008 0.04481 
 Residual             0.002574 0.05074 
Number of obs: 120, groups:  Genotype, 3

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)               0.16891    0.02748   2.38972   6.147 0.016262 *  
OriginTransplant         -0.02845    0.01310 114.00000  -2.172 0.031936 *  
SiteSS                   -0.05868    0.01310 114.00000  -4.479 1.79e-05 ***
OriginTransplant:SiteSS   0.06284    0.01853 114.00000   3.392 0.000954 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.238              
SiteSS      -0.238  0.500       
OrgnTrn:SSS  0.169 -0.707 -0.707

eta_squared(TLE_TP3.4.lme)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       9.01e-04 | [0.00, 1.00]
Site        |           0.07 | [0.01, 1.00]
Origin:Site |           0.09 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
TLE_TP3.4.res<-data.frame(summary(TLE_TP3.4.lme)$coefficients[-1,])
TLE_TP3.4.res$Predictor<-c("Origin", "Site", "Origin x Site")
TLE_TP3.4.res$EtaSq<-c(eta_squared(TLE_TP3.4.lme)$Eta2)
TLE_TP3.4.res$Response<-rep("Growth", nrow(TLE_TP3.4.res))
TLE_TP3.4.res$Timepoint<-rep("TP3_4", nrow(TLE_TP3.4.res))

Pairwise

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
TLE.pair_TP3.4<-emmeans(TLE_TP3.4.lme, pairwise~Origin | Site)
TLE.pair_TP3.4

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.169 0.0275 2.39  0.06737    0.270
 Transplant  0.140 0.0275 2.39  0.03892    0.242

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.110 0.0275 2.39  0.00869    0.212
 Transplant  0.145 0.0275 2.39  0.04308    0.246

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant   0.0285 0.0131 114   2.172  0.0319

Site = SS:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant  -0.0344 0.0131 114  -2.625  0.0098

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
TLE.pair.es_TP3.4<-data.frame(eff_size(TLE.pair_TP3.4, sigma=sigma(TLE_TP3.4.lme), edf=df.residual(TLE_TP3.4.lme)))

Since 'object' is a list, we are using the contrasts already present.

TLE.pair.es_TP3.4

TLE.pair.es_TP3.4<-TLE.pair.es_TP3.4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
TLE.pair_TP3.4.res<-merge(data.frame(TLE.pair_TP3.4$contrasts), TLE.pair.es_TP3.4[,-c(1)])
TLE.pair_TP3.4.res$Response<-rep("Growth", nrow(TLE.pair_TP3.4.res))
TLE.pair_TP3.4.res$TimeP<-rep("TP3_4", nrow(TLE.pair_TP3.4.res))

##Add Significance levels
TLE.pair_TP3.4.res$Sig<-ifelse(TLE.pair_TP3.4.res$p.value<0.001, "***",
                              ifelse(TLE.pair_TP3.4.res$p.value<0.01, "**",
                                     ifelse(TLE.pair_TP3.4.res$p.value<0.05, "*",
                                            ifelse(TLE.pair_TP3.4.res$p.value<0.1, "-", NA))))
TLE.pair_TP3.4.res

Plot Growth TP 3 to 4

##Summary statistics by Site and Origin
TP3.4_Growth.sum<-summarySE(Growth_TP3.4, measurevar="TLE_cm.day", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Growth Rate across Treatments
TP3.4_Growth.plot<-ggplot(TP3.4_Growth.sum, aes(x=Site, y=TLE_cm.day, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=TLE_cm.day-se, ymax=TLE_cm.day+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y=expression(paste('Growth Rate (cm day'^-1*")")), colour="Origin")+
  ylim(0, 0.25)+
  annotate("text", x=1, y=0.22, label="*", size=sig.sz, fontface="bold")+
  annotate("text", x=2, y=0.18, label="**", size=sig.sz, fontface="bold")+
  geom_bracket(xmin=1, xmax=2, y.position=0.01, label.size=levels.sz, label="***", inherit.aes = FALSE); TP3.4_Growth.plot

Percent Difference in Growth Rates

##Summary statistics by Genotype, Origin, and Site
TP3.4_Growth.sum<-summarySE(Growth_TP3.4, measurevar="TLE_cm.day", groupvars=c("Genotype", "Origin", "Site"), na.rm=TRUE)

##Calculate Percent Difference between Native and Transplant by Genotype
TP3.4_Growth.dif<-TP3.4_Growth.sum[which(TP3.4_Growth.sum$Origin=="Native"), c(1,3,5)]
TP3.4_Growth.dif<-TP3.4_Growth.dif %>% dplyr::rename(TLE_N = TLE_cm.day)
TP3.4_Growth.dif$TLE_T<-TP3.4_Growth.sum$TLE_cm.day[which(TP3.4_Growth.sum$Origin=="Transplant")]

##SS
TP3.4_Growth.dif_SS<-subset(TP3.4_Growth.dif, Site=="SS")
TP3.4_Growth.dif_SS$Perc.Inc<-((TP3.4_Growth.dif_SS$TLE_T-TP3.4_Growth.dif_SS$TLE_N)/TP3.4_Growth.dif_SS$TLE_N)*100
TP3.4_Growth.dif_SS

##Mean and Standard Error
mean(TP3.4_Growth.dif_SS$Perc.Inc)

[1] 34.22395

std.error(TP3.4_Growth.dif_SS$Perc.Inc)

[1] 18.61694

##KL
TP3.4_Growth.dif_KL<-subset(TP3.4_Growth.dif, Site=="KL")
TP3.4_Growth.dif_KL$Perc.Inc<-((TP3.4_Growth.dif_KL$TLE_N-TP3.4_Growth.dif_KL$TLE_T)/TP3.4_Growth.dif_KL$TLE_T)*100
TP3.4_Growth.dif_KL

##Mean and Standard Error
mean(TP3.4_Growth.dif_KL$Perc.Inc)

[1] 21.37099

std.error(TP3.4_Growth.dif_KL$Perc.Inc)

[1] 4.101183

Write Out Growth Data

##Growth Data dataframe
GrowthData<-Growth_TP1.2

##Specify Timepoints
GrowthData<-GrowthData %>% dplyr::rename(TP1.2_Ext_cm = Ext_cm, TP1.2_TLE_cm.day = TLE_cm.day)

##Merge with Timepoints 3 to 4
GrowthData<-merge(GrowthData, Growth_TP3.4, all=TRUE)

##Specify Timepoints
GrowthData<-GrowthData %>% dplyr::rename(TP3.4_Ext_cm = Ext_cm, TP3.4_TLE_cm.day = TLE_cm.day)

##Growth Data
write.csv(GrowthData, "Outputs/GrowthData.csv", row.names=FALSE)

Thermal Tolerance

Percent Retention

Subset Thermal Assay Data by Treatment

##Control
Therm_C<-subset(Thermal, Treat=="C")

##Heated
Therm_H<-subset(Thermal, Treat=="H")

Calculate Percent Retention

Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.

##Calculate averages for Control Treatment for each Site, Genotype, and Timepoint 
Therm_C<-na.omit(Therm_C)
names(Therm_C)

 [1] "ID"          "RandN"       "TimeP"       "Site"        "Genotype"    "Treat"      
 [7] "Treatment"   "Orig"        "Origin"      "Set"         "Site.Orig"   "SA_cm2"     
[13] "Chl_ug.cm2"  "Sym10.6_cm2" "Fv_Fm"      

Therm_C.a<-aggregate(Therm_C[,c(13:15)], list(Therm_C$Site, Therm_C$Genotype, Therm_C$TimeP, Therm_C$Origin), mean, na.action = na.omit)
names(Therm_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Origin")
names(Therm_C.a)[5:7]<-paste(names(Therm_C.a)[5:7], "C", sep="_")

##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, and Timepoint
Therm_H<-merge(Therm_H, Therm_C.a, all.x=TRUE )

##Calculate Proportion Retained for each Bleaching Metric relative to Control
Therm_H$Chl.prop<-round((Therm_H$Chl_ug.cm2/Therm_H$Chl_ug.cm2_C), 4)
Therm_H$Sym.prop<-round((Therm_H$Sym10.6_cm2/Therm_H$Sym10.6_cm2_C), 4)
Therm_H$PAM.prop<-round((Therm_H$Fv_Fm/Therm_H$Fv_Fm_C), 4)

##Set values >1 to 1
Therm_H$Chl.prop[which(Therm_H$Chl.prop>1)]<-1.0000
Therm_H$Sym.prop[which(Therm_H$Sym.prop>1)]<-1.0000
Therm_H$PAM.prop[which(Therm_H$PAM.prop>1)]<-1.0000

Initial Heat Assay

##Subset Initial Timepoint
Thermal_IN<-subset(Thermal, TimeP=="IN")
Thermal_IN$Site.Treat<-paste(Thermal_IN$Site, Thermal_IN$Treat, sep=".")
Thermal_IN$Site.Treat<-factor(Thermal_IN$Site.Treat, levels=c("KL.C", "KL.H", "SS.C", "SS.H"))


Therm_IN<-subset(Therm_H, TimeP=="IN")

Heated vs Control

Symbionts

#Check normality
hist(Thermal_IN$Sym10.6_cm2)

shapiro.test(Thermal_IN$Sym10.6_cm2)

    Shapiro-Wilk normality test

data:  Thermal_IN$Sym10.6_cm2
W = 0.90824, p-value = 0.001334

#Not normal

hist(log(Thermal_IN$Sym10.6_cm2+1))

shapiro.test(log(Thermal_IN$Sym10.6_cm2+1))

    Shapiro-Wilk normality test

data:  log(Thermal_IN$Sym10.6_cm2 + 1)
W = 0.95667, p-value = 0.07936

#Normal

##Model
#Function of Site and Treatment, with Genotype as a Random effect
Sym.lme_IN<-lmer(log(Sym10.6_cm2+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)

##Check residuals
Sym.lme_res_IN <- simulateResiduals(fittedModel = Sym.lme_IN, plot = F)
plot(Sym.lme_res_IN)

##Model results
summary(Sym.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym10.6_cm2 + 1) ~ Site * Treatment + (1 | Genotype)
   Data: Thermal_IN

REML criterion at convergence: -29.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.79325 -0.66013  0.08255  0.57806  2.12486 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0008222 0.02867 
 Residual             0.0230525 0.15183 
Number of obs: 47, groups:  Genotype, 3

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           0.38823    0.04685 14.06727   8.286 8.77e-07 ***
SiteSS                0.19508    0.06198 41.00308   3.147  0.00307 ** 
TreatmentHeat        -0.20370    0.06340 41.10985  -3.213  0.00256 ** 
SiteSS:TreatmentHeat -0.01724    0.08867 41.05799  -0.194  0.84677    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) SiteSS TrtmnH
SiteSS      -0.661              
TreatmentHt -0.647  0.489       
StSS:TrtmnH  0.462 -0.699 -0.715

eta_squared(Sym.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Site           |           0.30 | [0.12, 1.00]
Treatment      |           0.36 | [0.17, 1.00]
Site:Treatment |       9.20e-04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)

$emmeans
Treatment = Control:
 Site emmean     SE   df lower.CL upper.CL
 KL    0.388 0.0469 14.1   0.2878    0.489
 SS    0.583 0.0469 14.1   0.4829    0.684

Treatment = Heat:
 Site emmean     SE   df lower.CL upper.CL
 KL    0.185 0.0488 15.6   0.0808    0.288
 SS    0.362 0.0469 14.1   0.2619    0.463

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Treatment = Control:
 contrast estimate     SE   df t.ratio p.value
 KL - SS    -0.195 0.0620 41.0  -3.147  0.0031

Treatment = Heat:
 contrast estimate     SE   df t.ratio p.value
 KL - SS    -0.178 0.0635 41.1  -2.801  0.0077

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)

$emmeans
Site = KL:
 Treatment emmean     SE   df lower.CL upper.CL
 Control    0.388 0.0469 14.1   0.2878    0.489
 Heat       0.185 0.0488 15.6   0.0808    0.288

Site = SS:
 Treatment emmean     SE   df lower.CL upper.CL
 Control    0.583 0.0469 14.1   0.4829    0.684
 Heat       0.362 0.0469 14.1   0.2619    0.463

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast       estimate     SE   df t.ratio p.value
 Control - Heat    0.204 0.0635 41.1   3.209  0.0026

Site = SS:
 contrast       estimate     SE   df t.ratio p.value
 Control - Heat    0.221 0.0620 41.0   3.564  0.0009

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

##Save model results
Sym_IN.res<-data.frame(summary(Sym.lme_IN)$coefficients[-1,])

##Add pairwise
names(Sym_IN.res)<-names(data.frame(emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
Sym_IN.res<-rbind(Sym_IN.res,
data.frame(rbind(emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(Sym.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(Sym.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])

##Metadata
Sym_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
Sym_IN.res$Response<-rep("Symbionts", nrow(Sym_IN.res))
Sym_IN.res$EtaSq<-c(eta_squared(Sym.lme_IN)$Eta2, rep(NA, 4))

##Summary statistics by Site
IN_Sym.sum<-summarySE(Thermal_IN, measurevar="Sym10.6_cm2", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)

##Plot Average Symbionts across Treatments
IN_Sym.plot<-ggplot(IN_Sym.sum, aes(x=Site.Treat, y=Sym10.6_cm2, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Sym10.6_cm2-se, ymax=Sym10.6_cm2+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
  scale_shape_manual(values=c(1,16))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Treatment", y=expression(paste('Symbiont Density ('*10^6,'cells cm'^-2*")")), colour="Site")+
  ylim(0, 1)+
  annotate("text", x=.6, y=1, hjust=0, label="Site **", size=sig.sz, fontface="bold")+
  annotate("text", x=.6, y=.94, hjust=0, label="Treatment **", size=sig.sz, fontface="bold"); IN_Sym.plot

Chlorophyll

#Check normality
hist(Thermal_IN$Chl_ug.cm2)

shapiro.test(Thermal_IN$Chl_ug.cm2)

    Shapiro-Wilk normality test

data:  Thermal_IN$Chl_ug.cm2
W = 0.88059, p-value = 0.0002132

#Not normal

hist(log(Thermal_IN$Chl_ug.cm2+1))

shapiro.test(log(Thermal_IN$Chl_ug.cm2+1))

    Shapiro-Wilk normality test

data:  log(Thermal_IN$Chl_ug.cm2 + 1)
W = 0.92511, p-value = 0.005653

#Still not normal but less skewed

##Model
#Function of Site and Treatment, with Genotype as a Random effect
Chl.lme_IN<-lmer(log(Chl_ug.cm2+1)~Site*Treatment+(1|Genotype), data=Thermal_IN)

##Check residuals
Chl.lme_res_IN <- simulateResiduals(fittedModel = Chl.lme_IN, plot = F)
plot(Chl.lme_res_IN)

##Model results
summary(Chl.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl_ug.cm2 + 1) ~ Site * Treatment + (1 | Genotype)
   Data: Thermal_IN

REML criterion at convergence: -69.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.50589 -0.51517  0.07352  0.50169  2.35158 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.006137 0.07834 
 Residual             0.007930 0.08905 
Number of obs: 46, groups:  Genotype, 3

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           0.45648    0.05202  2.95779   8.774  0.00329 ** 
SiteSS                0.27810    0.03635 39.99255   7.650 2.38e-09 ***
TreatmentHeat        -0.37295    0.03721 40.00578 -10.024 1.80e-12 ***
SiteSS:TreatmentHeat -0.16117    0.05265 40.01278  -3.061  0.00393 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) SiteSS TrtmnH
SiteSS      -0.349              
TreatmentHt -0.341  0.489       
StSS:TrtmnH  0.241 -0.691 -0.707

eta_squared(Chl.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Site           |           0.58 | [0.41, 1.00]
Treatment      |           0.88 | [0.82, 1.00]
Site:Treatment |           0.19 | [0.04, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)

$emmeans
Treatment = Control:
 Site emmean     SE   df lower.CL upper.CL
 KL   0.4565 0.0520 2.97   0.2899    0.623
 SS   0.7346 0.0520 2.97   0.5680    0.901

Treatment = Heat:
 Site emmean     SE   df lower.CL upper.CL
 KL   0.0835 0.0526 3.11  -0.0808    0.248
 SS   0.2005 0.0526 3.11   0.0361    0.365

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Treatment = Control:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.278 0.0364 40  -7.650  <.0001

Treatment = Heat:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.117 0.0381 40  -3.070  0.0038

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)

$emmeans
Site = KL:
 Treatment emmean     SE   df lower.CL upper.CL
 Control   0.4565 0.0520 2.97   0.2899    0.623
 Heat      0.0835 0.0526 3.11  -0.0808    0.248

Site = SS:
 Treatment emmean     SE   df lower.CL upper.CL
 Control   0.7346 0.0520 2.97   0.5680    0.901
 Heat      0.2005 0.0526 3.11   0.0361    0.365

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast       estimate     SE df t.ratio p.value
 Control - Heat    0.373 0.0372 40  10.022  <.0001

Site = SS:
 contrast       estimate     SE df t.ratio p.value
 Control - Heat    0.534 0.0372 40  14.353  <.0001

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

##Save model results
Chl_IN.res<-data.frame(summary(Chl.lme_IN)$coefficients[-1,])

##Add pairwise
names(Chl_IN.res)<-names(data.frame(emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
Chl_IN.res<-rbind(Chl_IN.res,
data.frame(rbind(emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(Chl.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(Chl.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])

##Metadata
Chl_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
Chl_IN.res$Response<-rep("Chlorophyll", nrow(Chl_IN.res))
Chl_IN.res$EtaSq<-c(eta_squared(Chl.lme_IN)$Eta2, rep(NA, 4))

##Summary statistics by Site
IN_Chl.sum<-summarySE(Thermal_IN, measurevar="Chl_ug.cm2", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)

##Plot Average Chlorophyll across Treatments
IN_Chl.plot<-ggplot(IN_Chl.sum, aes(x=Site.Treat, y=Chl_ug.cm2, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Chl_ug.cm2-se, ymax=Chl_ug.cm2+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
  scale_shape_manual(values=c(1,16))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Treatment", y=expression(paste('Total Chlorophyll (\u03BCg cm'^-2*")")), colour="Site")+
  ylim(0, 1.25)+
  annotate("text", x=.6, y=1.2, hjust=0, label="Site ***", size=sig.sz, fontface="bold")+
  annotate("text", x=.6, y=1.12, hjust=0, label="Treatment ***", size=sig.sz, fontface="bold")+
  annotate("text", x=.6, y=1.04, hjust=0, label="Site x Treatment **", size=sig.sz, fontface="bold"); IN_Chl.plot

FvFm

#Check normality
hist(Thermal_IN$Fv_Fm)

shapiro.test(Thermal_IN$Fv_Fm)

    Shapiro-Wilk normality test

data:  Thermal_IN$Fv_Fm
W = 0.87006, p-value = 9.097e-05

#Not normal

hist(Thermal_IN$Fv_Fm^2)

shapiro.test(Thermal_IN$Fv_Fm^2)

    Shapiro-Wilk normality test

data:  Thermal_IN$Fv_Fm^2
W = 0.90968, p-value = 0.001489

#Still not normal but less skewed

##Model
#Function of Site and Treatment, with Genotype as a Random effect
PAM.lme_IN<-lmer(Fv_Fm^2~Site*Treatment+(1|Genotype), data=Thermal_IN)

boundary (singular) fit: see help('isSingular')

##Check residuals
PAM.lme_res_IN <- simulateResiduals(fittedModel = PAM.lme_IN, plot = F)
plot(PAM.lme_res_IN)

##Model results
summary(PAM.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Fv_Fm^2 ~ Site * Treatment + (1 | Genotype)
   Data: Thermal_IN

REML criterion at convergence: -156.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2454 -0.5025  0.1248  0.7046  1.5447 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 Genotype (Intercept) 8.057e-21 8.976e-11
 Residual             1.235e-03 3.514e-02
Number of obs: 47, groups:  Genotype, 3

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           0.38027    0.01015 43.00000  37.482  < 2e-16 ***
SiteSS                0.02899    0.01435 43.00000   2.021   0.0496 *  
TreatmentHeat        -0.07814    0.01467 43.00000  -5.327 3.45e-06 ***
SiteSS:TreatmentHeat  0.03376    0.02052 43.00000   1.645   0.1072    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) SiteSS TrtmnH
SiteSS      -0.707              
TreatmentHt -0.692  0.489       
StSS:TrtmnH  0.494 -0.699 -0.715
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(PAM.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Site           |           0.32 | [0.14, 1.00]
Treatment      |           0.45 | [0.27, 1.00]
Site:Treatment |           0.06 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)

$emmeans
Treatment = Control:
 Site emmean     SE   df lower.CL upper.CL
 KL    0.380 0.0101 21.6    0.359    0.401
 SS    0.409 0.0101 21.6    0.388    0.430

Treatment = Heat:
 Site emmean     SE   df lower.CL upper.CL
 KL    0.302 0.0106 23.5    0.280    0.324
 SS    0.365 0.0101 21.6    0.344    0.386

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Treatment = Control:
 contrast estimate     SE   df t.ratio p.value
 KL - SS   -0.0290 0.0143 41.0  -2.021  0.0499

Treatment = Heat:
 contrast estimate     SE   df t.ratio p.value
 KL - SS   -0.0627 0.0147 41.2  -4.268  0.0001

Degrees-of-freedom method: kenward-roger 

emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)

$emmeans
Site = KL:
 Treatment emmean     SE   df lower.CL upper.CL
 Control    0.380 0.0101 21.6    0.359    0.401
 Heat       0.302 0.0106 23.5    0.280    0.324

Site = SS:
 Treatment emmean     SE   df lower.CL upper.CL
 Control    0.409 0.0101 21.6    0.388    0.430
 Heat       0.365 0.0101 21.6    0.344    0.386

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast       estimate     SE   df t.ratio p.value
 Control - Heat   0.0781 0.0147 41.2   5.316  <.0001

Site = SS:
 contrast       estimate     SE   df t.ratio p.value
 Control - Heat   0.0444 0.0143 41.0   3.094  0.0036

Degrees-of-freedom method: kenward-roger 

##Save model results
PAM_IN.res<-data.frame(summary(PAM.lme_IN)$coefficients[-1,])

##Add pairwise
names(PAM_IN.res)<-names(data.frame(emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts))[-c(1:2)]
PAM_IN.res<-rbind(PAM_IN.res,
data.frame(rbind(emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts[1], emmeans(PAM.lme_IN, pairwise ~ Site | Treatment)$contrasts[2], emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)$contrasts[1], emmeans(PAM.lme_IN, pairwise ~ Treatment | Site)$contrasts[2]))[,-c(1:3)])

##Metadata
PAM_IN.res$Predictor<-c("Site", "Treatment", "Site x Treatment", "Control Site", "Heated Site", "KL Treatment", "SS Treatment")
PAM_IN.res$Response<-rep("FvFm", nrow(PAM_IN.res))
PAM_IN.res$EtaSq<-c(eta_squared(PAM.lme_IN)$Eta2, rep(NA, 4))

Not normal, but model residuals are OK.

##Summary statistics by Site
IN_PAM.sum<-summarySE(Thermal_IN, measurevar="Fv_Fm", groupvars=c("Site", "Treatment", "Site.Treat"), na.rm=TRUE)

##Plot Average FvFm across Treatments
IN_PAM.plot<-ggplot(IN_PAM.sum, aes(x=Site.Treat, y=Fv_Fm, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Fv_Fm-se, ymax=Fv_Fm+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Treatment), size=point.sz, position=position_dodge(width=0.5))+
  scale_shape_manual(values=c(1,16))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Treatment", y="Photochemical Efficiency (Fv/Fm)", colour="Site")+
  ylim(0.5, 0.655)+
  annotate("text", x=.6, y=.655, hjust=0, label="Site *", size=sig.sz, fontface="bold")+
  annotate("text", x=.6, y=.645, hjust=0, label="Treatment ***", size=sig.sz, fontface="bold"); IN_PAM.plot

Retention by Site

Symbionts

#Check normality
hist(Therm_IN$Sym.prop)

shapiro.test(Therm_IN$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_IN$Sym.prop
W = 0.93042, p-value = 0.1117

#Normal

##Model
#Function of Site, with Genotype as a Random effect
Tol_Sym.lme_IN<-lmer(Sym.prop~Site+(1|Genotype), data=Therm_IN)

##Check residuals
Tol_Sym.lme_res_IN <- simulateResiduals(fittedModel = Tol_Sym.lme_IN, plot = F)
plot(Tol_Sym.lme_res_IN)

##Model results
summary(Tol_Sym.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Site + (1 | Genotype)
   Data: Therm_IN

REML criterion at convergence: -2.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6802 -0.6578 -0.1521  0.6805  1.9852 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.07818  0.2796  
 Residual             0.03039  0.1743  
Number of obs: 23, groups:  Genotype, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)  0.46139    0.16983  2.21175   2.717   0.1014  
SiteSS       0.14316    0.07291 19.00328   1.964   0.0644 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
SiteSS -0.225

eta_squared(Tol_Sym.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Site      |           0.17 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_IN.res<-data.frame(summary(Tol_Sym.lme_IN)$coefficients)
names(Tol_Sym_IN.res)<-names(Sym_IN.res)[1:5]
Tol_Sym_IN.res$Predictor<-rep("Site", nrow(Tol_Sym_IN.res))
Tol_Sym_IN.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_IN.res))
Tol_Sym_IN.res$EtaSq<-c(eta_squared(Tol_Sym.lme_IN)$Eta2)

##Combine results
Sym_IN.res<-rbind(Sym_IN.res, Tol_Sym_IN.res[2,])

##Summary statistics by Site
IN_TolSym.sum<-summarySE(Therm_IN, measurevar="Sym.prop", groupvars=c("Site"), na.rm=TRUE)

##Plot Average Symbiont Retention across Treatments
IN_TolSym.plot<-ggplot(IN_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site", y="Symbiont Retention", colour="Site")+
  ylim(0, 1)+
  annotate("text", x=1.5, y=0.75, label="-", size=sig.sz, fontface="bold"); IN_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_IN$Chl.prop)

shapiro.test(Therm_IN$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_IN$Chl.prop
W = 0.98432, p-value = 0.9694

#Normal

##Model
#Function of Site, with Genotype as a Random effect
Tol_Chl.lme_IN<-lmer(Chl.prop~Site+(1|Genotype), data=Therm_IN)

##Check residuals
Tol_Chl.lme_res_IN <- simulateResiduals(fittedModel = Tol_Chl.lme_IN, plot = F)
plot(Tol_Chl.lme_res_IN)

##Model results
summary(Tol_Chl.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Chl.prop ~ Site + (1 | Genotype)
   Data: Therm_IN

REML criterion at convergence: -47.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3909 -0.5550 -0.1823  0.5284  1.9947 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.004576 0.06765 
 Residual             0.003436 0.05861 
Number of obs: 22, groups:  Genotype, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)  0.16552    0.04290  2.40867   3.858    0.045 *
SiteSS       0.04306    0.02514 18.06310   1.713    0.104  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
SiteSS -0.293

eta_squared(Tol_Chl.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Site      |           0.14 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_IN.res<-data.frame(summary(Tol_Chl.lme_IN)$coefficients)
names(Tol_Chl_IN.res)<-names(Chl_IN.res)[1:5]
Tol_Chl_IN.res$Predictor<-rep("Site", nrow(Tol_Chl_IN.res))
Tol_Chl_IN.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_IN.res))
Tol_Chl_IN.res$EtaSq<-c(eta_squared(Tol_Chl.lme_IN)$Eta2)

##Combine results
Chl_IN.res<-rbind(Chl_IN.res, Tol_Chl_IN.res[2,])

##Summary statistics by Site and Origin
IN_TolChl.sum<-summarySE(Therm_IN, measurevar="Chl.prop", groupvars=c("Site"), na.rm=TRUE)

##Plot Average Chlorophyll Retention across Treatments
IN_TolChl.plot<-ggplot(IN_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site", y="Chlorophyll Retention", colour="Site")+
  ylim(0, 0.35); IN_TolChl.plot

FvFm

#Check normality
hist(Therm_IN$PAM.prop)

shapiro.test(Therm_IN$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_IN$PAM.prop
W = 0.8289, p-value = 0.001154

#Not normal

hist(Therm_IN$PAM.prop^2)

shapiro.test(Therm_IN$PAM.prop^2)

    Shapiro-Wilk normality test

data:  Therm_IN$PAM.prop^2
W = 0.86491, p-value = 0.005132

##Still not normal

##Model
#Function of Site, with Genotype as a Random effect
Tol_PAM.lme_IN<-lmer(PAM.prop~Site+(1|Genotype), data=Therm_IN)

##Check residuals
Tol_PAM.lme_res_IN <- simulateResiduals(fittedModel = Tol_PAM.lme_IN, plot = F)
plot(Tol_PAM.lme_res_IN)

##Model results
summary(Tol_PAM.lme_IN)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Site + (1 | Genotype)
   Data: Therm_IN

REML criterion at convergence: -50.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7860 -0.2498  0.1225  0.5362  1.4591 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002221 0.04713 
 Residual             0.003538 0.05948 
Number of obs: 23, groups:  Genotype, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.89138    0.03262  2.77551  27.325 0.000185 ***
SiteSS       0.05230    0.02487 19.00103   2.103 0.049053 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
SiteSS -0.399

eta_squared(Tol_PAM.lme_IN)

# Effect Size for ANOVA (Type III)

Parameter | Eta2 (partial) |       95% CI
-----------------------------------------
Site      |           0.19 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_IN.res<-data.frame(summary(Tol_PAM.lme_IN)$coefficients)
names(Tol_PAM_IN.res)<-names(PAM_IN.res)[1:5]
Tol_PAM_IN.res$Predictor<-rep("Site", nrow(Tol_PAM_IN.res))
Tol_PAM_IN.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_IN.res))
Tol_PAM_IN.res$EtaSq<-c(eta_squared(Tol_PAM.lme_IN)$Eta2)

##Combine results
PAM_IN.res<-rbind(PAM_IN.res, Tol_PAM_IN.res[2,])

Not normal, but model residuals are OK.

##Summary statistics by Site and Origin
IN_TolPAM.sum<-summarySE(Therm_IN, measurevar="PAM.prop", groupvars=c("Site"), na.rm=TRUE)

##Plot Average FvFm Retention across Treatments
IN_TolPAM.plot<-ggplot(IN_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site", y="Photochemical Efficiency Retention", colour="Site")+
  ylim(0.6, 1)+
  annotate("text", x=1.5, y=1, label="*", size=sig.sz, fontface="bold"); IN_TolPAM.plot

Tolerance after Transplant

TP1

##Subset Timepoint 1
Therm_TP1<-subset(Therm_H, TimeP=="TP1")

Symbionts

#Check normality
hist(Therm_TP1$Sym.prop)

shapiro.test(Therm_TP1$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP1$Sym.prop
W = 0.91699, p-value = 0.002332

#Not normal

hist(log(Therm_TP1$Sym.prop+1))

shapiro.test(log(Therm_TP1$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP1$Sym.prop + 1)
W = 0.92258, p-value = 0.003671

##Still not normal 

##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP1<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP1)

##Check residuals
Tol_Sym.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP1, plot = F)
plot(Tol_Sym.lme_res_TP1)

##Model results
summary(Tol_Sym.lme_TP1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP1

REML criterion at convergence: -3.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2276 -0.5910 -0.0696  0.7547  1.8547 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.07291  0.2700  
 Residual             0.03714  0.1927  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)              0.35578    0.16552  2.38641   2.149   0.1436  
OriginTransplant         0.14490    0.07867 42.00000   1.842   0.0726 .
SiteSS                   0.06829    0.07867 42.00000   0.868   0.3903  
OriginTransplant:SiteSS -0.15161    0.11126 42.00000  -1.363   0.1803  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.238              
SiteSS      -0.238  0.500       
OrgnTrn:SSS  0.168 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP1)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |       4.34e-04 | [0.00, 1.00]
Origin:Site |           0.04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP1.res<-data.frame(summary(Tol_Sym.lme_TP1)$coefficients[-1,])
Tol_Sym_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP1.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP1)$Eta2)
Tol_Sym_TP1.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP1.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP1<-emmeans(Tol_Sym.lme_TP1, pairwise~Origin | Site)
Tol_Sym.pair_TP1

$emmeans
Site = KL:
 Origin     emmean    SE   df lower.CL upper.CL
 Native      0.356 0.166 2.39   -0.257    0.968
 Transplant  0.501 0.166 2.39   -0.112    1.113

Site = SS:
 Origin     emmean    SE   df lower.CL upper.CL
 Native      0.424 0.166 2.39   -0.188    1.036
 Transplant  0.417 0.166 2.39   -0.195    1.030

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant -0.14490 0.0787 42  -1.842  0.0726

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  0.00671 0.0787 42   0.085  0.9325

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP1<-data.frame(eff_size(Tol_Sym.pair_TP1, sigma=sigma(Tol_Sym.lme_TP1), edf=df.residual(Tol_Sym.lme_TP1)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Sym.pair.es_TP1

Tol_Sym.pair.es_TP1<-Tol_Sym.pair.es_TP1 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Sym.pair_TP1.res<-merge(data.frame(Tol_Sym.pair_TP1$contrasts), Tol_Sym.pair.es_TP1[,-c(1)])
Tol_Sym.pair_TP1.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP1.res))

##Add Significance levels
Tol_Sym.pair_TP1.res$Sig<-ifelse(Tol_Sym.pair_TP1.res$p.value<0.001, "***",
                              ifelse(Tol_Sym.pair_TP1.res$p.value<0.01, "**",
                                     ifelse(Tol_Sym.pair_TP1.res$p.value<0.05, "*",
                                            ifelse(Tol_Sym.pair_TP1.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP1.res

##Summary statistics by Site and Origin
TP1_TolSym.sum<-summarySE(Therm_TP1, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Symbiont Retention across Treatments
TP1_TolSym.plot<-ggplot(TP1_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1)+
  annotate("text", x=1, y=0.65, label="-", size=sig.sz, fontface="bold"); TP1_TolSym.plot

Chlorophyll

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site and Timepoint
Tol_Chl.pair<-emmeans(Tol_Chl.lme, pairwise~Origin | Site*TimeP)
Tol_Chl.pair

$emmeans
Site = KL, TimeP = TP1:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.1556 0.0408 41.3  0.07319    0.238
 Transplant 0.1697 0.0408 41.3  0.08727    0.252

Site = SS, TimeP = TP1:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.1849 0.0408 41.3  0.10249    0.267
 Transplant 0.1726 0.0408 41.3  0.09010    0.255

Site = KL, TimeP = TP2:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.3108 0.0408 41.3  0.22832    0.393
 Transplant 0.2723 0.0408 41.3  0.18988    0.355

Site = SS, TimeP = TP2:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.2456 0.0408 41.3  0.16318    0.328
 Transplant 0.2688 0.0408 41.3  0.18639    0.351

Site = KL, TimeP = TP3:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.5549 0.0408 41.3  0.47243    0.637
 Transplant 0.6116 0.0408 41.3  0.52910    0.694

Site = SS, TimeP = TP3:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.6275 0.0408 41.3  0.54508    0.710
 Transplant 0.6382 0.0408 41.3  0.55573    0.721

Site = KL, TimeP = TP4:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.1000 0.0408 41.3  0.01759    0.182
 Transplant 0.0858 0.0408 41.3  0.00330    0.168

Site = SS, TimeP = TP4:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.0723 0.0424 46.5 -0.01305    0.158
 Transplant 0.0895 0.0408 41.3  0.00707    0.172

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL, TimeP = TP1:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant  -0.0141 0.0533 173  -0.264  0.7919

Site = SS, TimeP = TP1:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant   0.0124 0.0533 173   0.233  0.8163

Site = KL, TimeP = TP2:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant   0.0384 0.0533 173   0.722  0.4715

Site = SS, TimeP = TP2:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant  -0.0232 0.0533 173  -0.436  0.6635

Site = KL, TimeP = TP3:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant  -0.0567 0.0533 173  -1.064  0.2889

Site = SS, TimeP = TP3:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant  -0.0107 0.0533 173  -0.200  0.8417

Site = KL, TimeP = TP4:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant   0.0143 0.0533 173   0.268  0.7888

Site = SS, TimeP = TP4:
 contrast            estimate     SE  df t.ratio p.value
 Native - Transplant  -0.0172 0.0545 173  -0.316  0.7523

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es<-data.frame(eff_size(Tol_Chl.pair, sigma=sigma(Tol_Chl.lme), edf=df.residual(Tol_Chl.lme)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Chl.pair.es

Tol_Chl.pair.es<-Tol_Chl.pair.es %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Chl.pair.res<-merge(data.frame(Tol_Chl.pair$contrasts), Tol_Chl.pair.es[,-c(1)])
Tol_Chl.pair.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP1<-emmeans(Tol_Chl.lme_TP1, pairwise~Origin | Site)
Tol_Chl.pair_TP1

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.156 0.0657 2.27  -0.0970    0.408
 Transplant  0.170 0.0657 2.27  -0.0829    0.422

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.185 0.0657 2.27  -0.0677    0.438
 Transplant  0.173 0.0657 2.27  -0.0801    0.425

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  -0.0141 0.0266 42  -0.529  0.5997

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant   0.0124 0.0266 42   0.466  0.6439

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP1<-data.frame(eff_size(Tol_Chl.pair_TP1, sigma=sigma(Tol_Chl.lme_TP1), edf=df.residual(Tol_Chl.lme_TP1)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Chl.pair.es_TP1

Tol_Chl.pair.es_TP1<-Tol_Chl.pair.es_TP1 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Chl.pair_TP1.res<-merge(data.frame(Tol_Chl.pair_TP1$contrasts), Tol_Chl.pair.es_TP1[,-c(1)])
Tol_Chl.pair_TP1.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP1.res))

##Add Significance levels
Tol_Chl.pair_TP1.res$Sig<-ifelse(Tol_Chl.pair_TP1.res$p.value<0.001, "***",
                              ifelse(Tol_Chl.pair_TP1.res$p.value<0.01, "**",
                                     ifelse(Tol_Chl.pair_TP1.res$p.value<0.05, "*",
                                            ifelse(Tol_Chl.pair_TP1.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP1.res

##Summary statistics by Site and Origin
TP1_TolChl.sum<-summarySE(Therm_TP1, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Chlorophyll Retention across Treatments
TP1_TolChl.plot<-ggplot(TP1_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 0.35); TP1_TolChl.plot

FvFm

#Check normality
hist(Therm_TP1$PAM.prop)

shapiro.test(Therm_TP1$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP1$PAM.prop
W = 0.84073, p-value = 1.245e-05

#Not normal

hist(Therm_TP1$PAM.prop^2)

shapiro.test(Therm_TP1$PAM.prop^2)

    Shapiro-Wilk normality test

data:  Therm_TP1$PAM.prop^2
W = 0.88155, p-value = 0.0001681

##Still not normal

##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP1<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP1)

##Check residuals
Tol_PAM.lme_res_TP1 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP1, plot = F)
plot(Tol_PAM.lme_res_TP1)

##Model results
summary(Tol_PAM.lme_TP1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP1

REML criterion at convergence: -84.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3058 -0.3498  0.1733  0.6717  1.7818 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.002735 0.05230 
 Residual             0.006214 0.07883 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.86181    0.03781  3.74555  22.794 3.72e-05 ***
OriginTransplant         0.06813    0.03218 42.00000   2.117   0.0402 *  
SiteSS                  -0.01277    0.03218 42.00000  -0.397   0.6936    
OriginTransplant:SiteSS -0.07942    0.04551 42.00000  -1.745   0.0883 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.426              
SiteSS      -0.426  0.500       
OrgnTrn:SSS  0.301 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP1)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |           0.11 | [0.01, 1.00]
Origin:Site |           0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP1.res<-data.frame(summary(Tol_PAM.lme_TP1)$coefficients[-1,])
Tol_PAM_TP1.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP1.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP1)$Eta2)
Tol_PAM_TP1.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP1.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP1<-emmeans(Tol_PAM.lme_TP1, pairwise~Origin | Site)
Tol_PAM.pair_TP1

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.862 0.0378 3.75    0.754    0.970
 Transplant  0.930 0.0378 3.75    0.822    1.038

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.849 0.0378 3.75    0.741    0.957
 Transplant  0.838 0.0378 3.75    0.730    0.946

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  -0.0681 0.0322 42  -2.117  0.0402

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant   0.0113 0.0322 42   0.351  0.7276

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP1<-data.frame(eff_size(Tol_PAM.pair_TP1, sigma=sigma(Tol_PAM.lme_TP1), edf=df.residual(Tol_PAM.lme_TP1)))

Since 'object' is a list, we are using the contrasts already present.

Tol_PAM.pair.es_TP1

Tol_PAM.pair.es_TP1<-Tol_PAM.pair.es_TP1 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_PAM.pair_TP1.res<-merge(data.frame(Tol_PAM.pair_TP1$contrasts), Tol_PAM.pair.es_TP1[,-c(1)])
Tol_PAM.pair_TP1.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP1.res))

##Add Significance levels
Tol_PAM.pair_TP1.res$Sig<-ifelse(Tol_PAM.pair_TP1.res$p.value<0.001, "***",
                              ifelse(Tol_PAM.pair_TP1.res$p.value<0.01, "**",
                                     ifelse(Tol_PAM.pair_TP1.res$p.value<0.05, "*",
                                            ifelse(Tol_PAM.pair_TP1.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP1.res

##Summary statistics by Site and Origin
TP1_TolPAM.sum<-summarySE(Therm_TP1, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average FvFm Retention across Treatments
TP1_TolPAM.plot<-ggplot(TP1_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.6, 1)+
  annotate("text", x=1, y=.96, label="*", size=sig.sz, fontface="bold"); TP1_TolPAM.plot

Save Results

##Combine Results
Tol_TP1.res<-rbind(Tol_PAM_TP1.res, Tol_Sym_TP1.res, Tol_Chl_TP1.res) 
Tol.pair.res_TP1<-rbind(Tol_PAM.pair_TP1.res, Tol_Sym.pair_TP1.res, Tol_Chl.pair_TP1.res)

##Add Timepoint
Tol_TP1.res$TimeP<-rep("TP1", nrow(Tol_TP1.res))
Tol.pair.res_TP1$TimeP<-rep("TP1", nrow(Tol.pair.res_TP1))

TP2

##Subset Timepoint 2
Therm_TP2<-subset(Therm_H, TimeP=="TP2")

Symbionts

#Check normality
hist(Therm_TP2$Sym.prop)

shapiro.test(Therm_TP2$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP2$Sym.prop
W = 0.89527, p-value = 0.0004439

#Not normal

hist(Therm_TP2$Sym.prop^2)

shapiro.test(Therm_TP2$Sym.prop^2)

    Shapiro-Wilk normality test

data:  Therm_TP2$Sym.prop^2
W = 0.8757, p-value = 0.0001129

##Still not normal

##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP2<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP2)

##Check residuals
Tol_Sym.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP2, plot = F)
plot(Tol_Sym.lme_res_TP2)

##Model results
summary(Tol_Sym.lme_TP2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP2

REML criterion at convergence: -3.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.20381 -0.54227  0.01202  0.47665  2.09092 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.04681  0.2164  
 Residual             0.03767  0.1941  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)              0.805833   0.136905  2.613359   5.886   0.0142 *
OriginTransplant        -0.107900   0.079237 42.000000  -1.362   0.1805  
SiteSS                  -0.178425   0.079237 42.000000  -2.252   0.0296 *
OriginTransplant:SiteSS  0.005183   0.112058 42.000000   0.046   0.9633  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.289              
SiteSS      -0.289  0.500       
OrgnTrn:SSS  0.205 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP2)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.08 | [0.00, 1.00]
Site        |           0.19 | [0.04, 1.00]
Origin:Site |       5.09e-05 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP2.res<-data.frame(summary(Tol_Sym.lme_TP2)$coefficients[-1,])
Tol_Sym_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP2.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP2)$Eta2)
Tol_Sym_TP2.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP2.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP2<-emmeans(Tol_Sym.lme_TP2, pairwise~Origin | Site)
Tol_Sym.pair_TP2

$emmeans
Site = KL:
 Origin     emmean    SE   df lower.CL upper.CL
 Native      0.806 0.137 2.61   0.3313    1.280
 Transplant  0.698 0.137 2.61   0.2234    1.172

Site = SS:
 Origin     emmean    SE   df lower.CL upper.CL
 Native      0.627 0.137 2.61   0.1529    1.102
 Transplant  0.525 0.137 2.61   0.0502    0.999

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant    0.108 0.0792 42   1.362  0.1805

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant    0.103 0.0792 42   1.296  0.2019

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP2<-data.frame(eff_size(Tol_Sym.pair_TP2, sigma=sigma(Tol_Sym.lme_TP2), edf=df.residual(Tol_Sym.lme_TP2)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Sym.pair.es_TP2

Tol_Sym.pair.es_TP2<-Tol_Sym.pair.es_TP2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Sym.pair_TP2.res<-merge(data.frame(Tol_Sym.pair_TP2$contrasts), Tol_Sym.pair.es_TP2[,-c(1)])
Tol_Sym.pair_TP2.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP2.res))

##Add Significance levels
Tol_Sym.pair_TP2.res$Sig<-ifelse(Tol_Sym.pair_TP2.res$p.value<0.001, "***",
                              ifelse(Tol_Sym.pair_TP2.res$p.value<0.01, "**",
                                     ifelse(Tol_Sym.pair_TP2.res$p.value<0.05, "*",
                                            ifelse(Tol_Sym.pair_TP2.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP2.res

##Summary statistics by Site and Origin
TP2_TolSym.sum<-summarySE(Therm_TP2, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP2_TolSym.plot<-ggplot(TP2_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1)+
 geom_bracket(xmin=1, xmax=2, y.position=.04, label.size=levels.sz, label="*", inherit.aes = FALSE); TP2_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP2$Chl.prop)

shapiro.test(Therm_TP2$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP2$Chl.prop
W = 0.97677, p-value = 0.4523

#Normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP2<-lmer(Chl.prop~Origin*Site+(1|Genotype), data=Therm_TP2)

##Check residuals
Tol_Chl.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP2, plot = F)
plot(Tol_Chl.lme_res_TP2)

##Model results
summary(Tol_Chl.lme_TP2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Chl.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP2

REML criterion at convergence: -81

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0024 -0.7622 -0.1174  0.6824  2.9666 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.006146 0.07840 
 Residual             0.006538 0.08086 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)  
(Intercept)              0.31077    0.05093  2.81331   6.102   0.0106 *
OriginTransplant        -0.03844    0.03301 42.00000  -1.165   0.2508  
SiteSS                  -0.06514    0.03301 42.00000  -1.973   0.0550 .
OriginTransplant:SiteSS  0.06166    0.04668 42.00000   1.321   0.1937  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.324              
SiteSS      -0.324  0.500       
OrgnTrn:SSS  0.229 -0.707 -0.707

eta_squared(Tol_Chl.lme_TP2)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       2.53e-03 | [0.00, 1.00]
Site        |           0.05 | [0.00, 1.00]
Origin:Site |           0.04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP2.res<-data.frame(summary(Tol_Chl.lme_TP2)$coefficients[-1,])
Tol_Chl_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP2.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP2)$Eta2)
Tol_Chl_TP2.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP2.res))

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP2<-emmeans(Tol_Chl.lme_TP2, pairwise~Origin | Site)
Tol_Chl.pair_TP2

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.311 0.0509 2.81   0.1424    0.479
 Transplant  0.272 0.0509 2.81   0.1040    0.441

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.246 0.0509 2.81   0.0773    0.414
 Transplant  0.269 0.0509 2.81   0.1005    0.437

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate    SE df t.ratio p.value
 Native - Transplant   0.0384 0.033 42   1.165  0.2508

Site = SS:
 contrast            estimate    SE df t.ratio p.value
 Native - Transplant  -0.0232 0.033 42  -0.703  0.4857

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP2<-data.frame(eff_size(Tol_Chl.pair_TP2, sigma=sigma(Tol_Chl.lme_TP2), edf=df.residual(Tol_Chl.lme_TP2)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Chl.pair.es_TP2

Tol_Chl.pair.es_TP2<-Tol_Chl.pair.es_TP2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Chl.pair_TP2.res<-merge(data.frame(Tol_Chl.pair_TP2$contrasts), Tol_Chl.pair.es_TP2[,-c(1)])
Tol_Chl.pair_TP2.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP2.res))

##Add Significance levels
Tol_Chl.pair_TP2.res$Sig<-ifelse(Tol_Chl.pair_TP2.res$p.value<0.001, "***",
                              ifelse(Tol_Chl.pair_TP2.res$p.value<0.01, "**",
                                     ifelse(Tol_Chl.pair_TP2.res$p.value<0.05, "*",
                                            ifelse(Tol_Chl.pair_TP2.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP2.res

##Summary statistics by Site and Origin
TP2_TolChl.sum<-summarySE(Therm_TP2, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP2_TolChl.plot<-ggplot(TP2_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 0.35)+
  geom_bracket(xmin=1, xmax=2, y.position=.02, label.size=levels.sz, label="-", inherit.aes = FALSE); TP2_TolChl.plot

FvFm

#Check normality
hist(Therm_TP2$PAM.prop)

shapiro.test(Therm_TP2$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP2$PAM.prop
W = 0.73527, p-value = 5.896e-08

#Not normal

hist(Therm_TP2$PAM.prop^2)

shapiro.test(Therm_TP2$PAM.prop^2)

    Shapiro-Wilk normality test

data:  Therm_TP2$PAM.prop^2
W = 0.80435, p-value = 1.633e-06

##Still not normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP2<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP2)

##Check residuals
Tol_PAM.lme_res_TP2 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP2, plot = F)
plot(Tol_PAM.lme_res_TP2)

##Model results
summary(Tol_PAM.lme_TP2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP2

REML criterion at convergence: -102.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8326 -0.3753 -0.0042  0.5350  1.7496 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.001675 0.04093 
 Residual             0.004192 0.06475 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.85622    0.03013  3.92146  28.419 1.09e-05 ***
OriginTransplant         0.07526    0.02643 42.00000   2.847  0.00680 ** 
SiteSS                   0.06487    0.02643 42.00000   2.454  0.01834 *  
OriginTransplant:SiteSS -0.10206    0.03738 42.00000  -2.730  0.00921 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.439              
SiteSS      -0.439  0.500       
OrgnTrn:SSS  0.310 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP2)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |           0.01 | [0.00, 1.00]
Origin:Site |           0.15 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP2.res<-data.frame(summary(Tol_PAM.lme_TP2)$coefficients[-1,])
Tol_PAM_TP2.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP2.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP2)$Eta2)
Tol_PAM_TP2.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP2.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP2<-emmeans(Tol_PAM.lme_TP2, pairwise~Origin | Site)
Tol_PAM.pair_TP2

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.856 0.0301 3.92    0.772    0.941
 Transplant  0.931 0.0301 3.92    0.847    1.016

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.921 0.0301 3.92    0.837    1.005
 Transplant  0.894 0.0301 3.92    0.810    0.979

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  -0.0753 0.0264 42  -2.847  0.0068

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant   0.0268 0.0264 42   1.014  0.3164

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP2<-data.frame(eff_size(Tol_PAM.pair_TP2, sigma=sigma(Tol_PAM.lme_TP2), edf=df.residual(Tol_PAM.lme_TP2)))

Since 'object' is a list, we are using the contrasts already present.

Tol_PAM.pair.es_TP2

Tol_PAM.pair.es_TP2<-Tol_PAM.pair.es_TP2 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_PAM.pair_TP2.res<-merge(data.frame(Tol_PAM.pair_TP2$contrasts), Tol_PAM.pair.es_TP2[,-c(1)])
Tol_PAM.pair_TP2.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP2.res))

##Add Significance levels
Tol_PAM.pair_TP2.res$Sig<-ifelse(Tol_PAM.pair_TP2.res$p.value<0.001, "***",
                              ifelse(Tol_PAM.pair_TP2.res$p.value<0.01, "**",
                                     ifelse(Tol_PAM.pair_TP2.res$p.value<0.05, "*",
                                            ifelse(Tol_PAM.pair_TP2.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP2.res

Save Results

##Combine Results
Tol_TP2.res<-rbind(Tol_PAM_TP2.res, Tol_Sym_TP2.res, Tol_Chl_TP2.res)
Tol.pair.res_TP2<-rbind(Tol_PAM.pair_TP2.res, Tol_Sym.pair_TP2.res, Tol_Chl.pair_TP2.res)

##Add Timepoint
Tol_TP2.res$TimeP<-rep("TP2", nrow(Tol_TP2.res))
Tol.pair.res_TP2$TimeP<-rep("TP2", nrow(Tol.pair.res_TP2))

TP3

##Subset Timepoint 3
Therm_TP3<-subset(Therm_H, TimeP=="TP3")

Symbionts

#Check normality
hist(Therm_TP3$Sym.prop)

shapiro.test(Therm_TP3$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP3$Sym.prop
W = 0.78653, p-value = 6.531e-07

#Not normal

hist(Therm_TP3$Sym.prop^2)

shapiro.test(Therm_TP3$Sym.prop^2)

    Shapiro-Wilk normality test

data:  Therm_TP3$Sym.prop^2
W = 0.78554, p-value = 6.213e-07

##Still not normal but less skewed

##Model with no transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP3<-lmer(Sym.prop~Origin*Site+(1|Genotype), data=Therm_TP3)

##Check residuals
Tol_Sym.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP3, plot = F)
plot(Tol_Sym.lme_res_TP3)

##Model results
summary(Tol_Sym.lme_TP3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Sym.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP3

REML criterion at convergence: -54.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1418 -0.5595  0.3785  0.6079  1.6260 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 0.0003012 0.01735 
 Residual             0.0132967 0.11531 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.82201    0.03476 16.65879  23.646 2.98e-14 ***
OriginTransplant         0.07462    0.04708 42.00000   1.585   0.1204    
SiteSS                   0.10814    0.04708 42.00000   2.297   0.0267 *  
OriginTransplant:SiteSS -0.05670    0.06658 42.00000  -0.852   0.3992    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.677              
SiteSS      -0.677  0.500       
OrgnTrn:SSS  0.479 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP3)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.04 | [0.00, 1.00]
Site        |           0.12 | [0.01, 1.00]
Origin:Site |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP3.res<-data.frame(summary(Tol_Sym.lme_TP3)$coefficients[-1,])
Tol_Sym_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP3.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP3)$Eta2)
Tol_Sym_TP3.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP3.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP3<-emmeans(Tol_Sym.lme_TP3, pairwise~Origin | Site)
Tol_Sym.pair_TP3

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.822 0.0348 16.7    0.749    0.895
 Transplant  0.897 0.0348 16.7    0.823    0.970

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.930 0.0348 16.7    0.857    1.004
 Transplant  0.948 0.0348 16.7    0.875    1.022

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  -0.0746 0.0471 42  -1.585  0.1204

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  -0.0179 0.0471 42  -0.381  0.7053

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP3<-data.frame(eff_size(Tol_Sym.pair_TP3, sigma=sigma(Tol_Sym.lme_TP3), edf=df.residual(Tol_Sym.lme_TP3)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Sym.pair.es_TP3

Tol_Sym.pair.es_TP3<-Tol_Sym.pair.es_TP3 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Sym.pair_TP3.res<-merge(data.frame(Tol_Sym.pair_TP3$contrasts), Tol_Sym.pair.es_TP3[,-c(1)])
Tol_Sym.pair_TP3.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP3.res))

##Add Significance levels
Tol_Sym.pair_TP3.res$Sig<-ifelse(Tol_Sym.pair_TP3.res$p.value<0.001, "***",
                              ifelse(Tol_Sym.pair_TP3.res$p.value<0.01, "**",
                                     ifelse(Tol_Sym.pair_TP3.res$p.value<0.05, "*",
                                            ifelse(Tol_Sym.pair_TP3.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP3.res

##Summary statistics by Site and Origin
TP3_TolSym.sum<-summarySE(Therm_TP3, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP3_TolSym.plot<-ggplot(TP3_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1)+
 geom_bracket(xmin=1, xmax=2, y.position=.04, label.size=levels.sz, label="*", inherit.aes = FALSE); TP3_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP3$Chl.prop)

shapiro.test(Therm_TP3$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP3$Chl.prop
W = 0.89534, p-value = 0.0004461

#Not normal

hist(log(Therm_TP3$Chl.prop+1))

shapiro.test(log(Therm_TP3$Chl.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP3$Chl.prop + 1)
W = 0.9211, p-value = 0.003251

##Still not normal but less skewed

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP3<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP3)

boundary (singular) fit: see help('isSingular')

##Check residuals
Tol_Chl.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP3, plot = F)
plot(Tol_Chl.lme_res_TP3)

##Model results
summary(Tol_Chl.lme_TP3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP3

REML criterion at convergence: -47.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7095 -0.6080 -0.1405  0.4999  1.7650 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.00000  0.0000  
 Residual             0.01594  0.1262  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.43851    0.03644 44.00000  12.032 1.65e-15 ***
OriginTransplant         0.03181    0.05154 44.00000   0.617    0.540    
SiteSS                   0.03852    0.05154 44.00000   0.747    0.459    
OriginTransplant:SiteSS -0.02511    0.07289 44.00000  -0.344    0.732    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.707              
SiteSS      -0.707  0.500       
OrgnTrn:SSS  0.500 -0.707 -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP3)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       6.31e-03 | [0.00, 1.00]
Site        |           0.01 | [0.00, 1.00]
Origin:Site |       2.69e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP3.res<-data.frame(summary(Tol_Chl.lme_TP3)$coefficients[-1,])
Tol_Chl_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP3.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP3)$Eta2)
Tol_Chl_TP3.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP3.res))

Not normal, but model residuals are OK.

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP3<-emmeans(Tol_Chl.lme_TP3, pairwise~Origin | Site)
Tol_Chl.pair_TP3

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.439 0.0364 22.4    0.363    0.514
 Transplant  0.470 0.0364 22.4    0.395    0.546

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.477 0.0364 22.4    0.402    0.553
 Transplant  0.484 0.0364 22.4    0.408    0.559

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant -0.03181 0.0515 42  -0.617  0.5404

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant -0.00671 0.0515 42  -0.130  0.8971

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP3<-data.frame(eff_size(Tol_Chl.pair_TP3, sigma=sigma(Tol_Chl.lme_TP3), edf=df.residual(Tol_Chl.lme_TP3)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Chl.pair.es_TP3

Tol_Chl.pair.es_TP3<-Tol_Chl.pair.es_TP3 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Chl.pair_TP3.res<-merge(data.frame(Tol_Chl.pair_TP3$contrasts), Tol_Chl.pair.es_TP3[,-c(1)])
Tol_Chl.pair_TP3.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP3.res))

##Add Significance levels
Tol_Chl.pair_TP3.res$Sig<-ifelse(Tol_Chl.pair_TP3.res$p.value<0.001, "***",
                              ifelse(Tol_Chl.pair_TP3.res$p.value<0.01, "**",
                                     ifelse(Tol_Chl.pair_TP3.res$p.value<0.05, "*",
                                            ifelse(Tol_Chl.pair_TP3.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP3.res

##Summary statistics by Site and Origin
TP3_TolChl.sum<-summarySE(Therm_TP3, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP3_TolChl.plot<-ggplot(TP3_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 1); TP3_TolChl.plot

FvFm

#Check normality
hist(Therm_TP3$PAM.prop)

shapiro.test(Therm_TP3$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP3$PAM.prop
W = 0.96358, p-value = 0.1408

#Normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP3<-lmer(PAM.prop~Origin*Site+(1|Genotype), data=Therm_TP3)

##Check residuals
Tol_PAM.lme_res_TP3 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP3, plot = F)
plot(Tol_PAM.lme_res_TP3)

##Model results
summary(Tol_PAM.lme_TP3)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP3

REML criterion at convergence: -172.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.29399 -0.71463  0.01706  0.83761  1.77082 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Genotype (Intercept) 3.177e-05 0.005637
 Residual             9.061e-04 0.030101
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                         Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)              0.954325   0.009279 14.520200 102.850   <2e-16 ***
OriginTransplant         0.009025   0.012289 42.000000   0.734    0.467    
SiteSS                   0.003975   0.012289 42.000000   0.323    0.748    
OriginTransplant:SiteSS -0.019133   0.017379 42.000000  -1.101    0.277    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.662              
SiteSS      -0.662  0.500       
OrgnTrn:SSS  0.468 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP3)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       9.25e-05 | [0.00, 1.00]
Site        |       9.76e-03 | [0.00, 1.00]
Origin:Site |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP3.res<-data.frame(summary(Tol_PAM.lme_TP3)$coefficients[-1,])
Tol_PAM_TP3.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP3.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP3)$Eta2)
Tol_PAM_TP3.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP3.res))

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP3<-emmeans(Tol_PAM.lme_TP3, pairwise~Origin | Site)
Tol_PAM.pair_TP3

$emmeans
Site = KL:
 Origin     emmean      SE   df lower.CL upper.CL
 Native      0.954 0.00928 14.5    0.934    0.974
 Transplant  0.963 0.00928 14.5    0.944    0.983

Site = SS:
 Origin     emmean      SE   df lower.CL upper.CL
 Native      0.958 0.00928 14.5    0.938    0.978
 Transplant  0.948 0.00928 14.5    0.928    0.968

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant -0.00902 0.0123 42  -0.734  0.4668

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  0.01011 0.0123 42   0.823  0.4154

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP3<-data.frame(eff_size(Tol_PAM.pair_TP3, sigma=sigma(Tol_PAM.lme_TP3), edf=df.residual(Tol_PAM.lme_TP3)))

Since 'object' is a list, we are using the contrasts already present.

Tol_PAM.pair.es_TP3

Tol_PAM.pair.es_TP3<-Tol_PAM.pair.es_TP3 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_PAM.pair_TP3.res<-merge(data.frame(Tol_PAM.pair_TP3$contrasts), Tol_PAM.pair.es_TP3[,-c(1)])
Tol_PAM.pair_TP3.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP3.res))

##Add Significance levels
Tol_PAM.pair_TP3.res$Sig<-ifelse(Tol_PAM.pair_TP3.res$p.value<0.001, "***",
                              ifelse(Tol_PAM.pair_TP3.res$p.value<0.01, "**",
                                     ifelse(Tol_PAM.pair_TP3.res$p.value<0.05, "*",
                                            ifelse(Tol_PAM.pair_TP3.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP3.res

##Summary statistics by Site and Origin
TP3_TolPAM.sum<-summarySE(Therm_TP3, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP3_TolPAM.plot<-ggplot(TP3_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.7, 1); TP3_TolPAM.plot

Save Results

##Combine Results
Tol_TP3.res<-rbind(Tol_PAM_TP3.res, Tol_Sym_TP3.res, Tol_Chl_TP3.res)
Tol.pair.res_TP3<-rbind(Tol_PAM.pair_TP3.res, Tol_Sym.pair_TP3.res, Tol_Chl.pair_TP3.res)

##Add Timepoint
Tol_TP3.res$TimeP<-rep("TP3", nrow(Tol_TP3.res))
Tol.pair.res_TP3$TimeP<-rep("TP3", nrow(Tol.pair.res_TP3))

TP4

##Subset Timepoint 4
Therm_TP4<-subset(Therm_H, TimeP=="TP4")

Symbionts

#Check normality
hist(Therm_TP4$Sym.prop)

shapiro.test(Therm_TP4$Sym.prop)

    Shapiro-Wilk normality test

data:  Therm_TP4$Sym.prop
W = 0.93453, p-value = 0.01005

#Not normal

hist(log(Therm_TP4$Sym.prop+1))

shapiro.test(log(Therm_TP4$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP4$Sym.prop + 1)
W = 0.95458, p-value = 0.061

#Normal

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Sym.lme_TP4<-lmer(log(Sym.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP4)

##Check residuals
Tol_Sym.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_Sym.lme_TP4, plot = F)
plot(Tol_Sym.lme_res_TP4)

##Model results
summary(Tol_Sym.lme_TP4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Sym.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP4

REML criterion at convergence: -44.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5412 -0.7855 -0.1534  0.5518  2.2412 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.001295 0.03598 
 Residual             0.016289 0.12763 
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.24607    0.04230  9.94436   5.818 0.000173 ***
OriginTransplant         0.04098    0.05210 42.00002   0.787 0.435975    
SiteSS                  -0.03092    0.05210 42.00002  -0.593 0.556032    
OriginTransplant:SiteSS -0.03256    0.07369 42.00002  -0.442 0.660816    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.616              
SiteSS      -0.616  0.500       
OrgnTrn:SSS  0.436 -0.707 -0.707

eta_squared(Tol_Sym.lme_TP4)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |           0.01 | [0.00, 1.00]
Site        |           0.04 | [0.00, 1.00]
Origin:Site |       4.63e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_TP4.res<-data.frame(summary(Tol_Sym.lme_TP4)$coefficients[-1,])
Tol_Sym_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Sym_TP4.res$EtaSq<-c(eta_squared(Tol_Sym.lme_TP4)$Eta2)
Tol_Sym_TP4.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym_TP4.res))

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Sym.pair_TP4<-emmeans(Tol_Sym.lme_TP4, pairwise~Origin | Site)
Tol_Sym.pair_TP4

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.246 0.0423 9.94    0.152    0.340
 Transplant  0.287 0.0423 9.94    0.193    0.381

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.215 0.0423 9.94    0.121    0.309
 Transplant  0.224 0.0423 9.94    0.129    0.318

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant -0.04098 0.0521 42  -0.787  0.4360

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant -0.00842 0.0521 42  -0.162  0.8724

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Sym.pair.es_TP4<-data.frame(eff_size(Tol_Sym.pair_TP4, sigma=sigma(Tol_Sym.lme_TP4), edf=df.residual(Tol_Sym.lme_TP4)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Sym.pair.es_TP4

Tol_Sym.pair.es_TP4<-Tol_Sym.pair.es_TP4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Sym.pair_TP4.res<-merge(data.frame(Tol_Sym.pair_TP4$contrasts), Tol_Sym.pair.es_TP4[,-c(1)])
Tol_Sym.pair_TP4.res$Response<-rep("Symbiont Retention", nrow(Tol_Sym.pair_TP4.res))

##Add Significance levels
Tol_Sym.pair_TP4.res$Sig<-ifelse(Tol_Sym.pair_TP4.res$p.value<0.001, "***",
                              ifelse(Tol_Sym.pair_TP4.res$p.value<0.01, "**",
                                     ifelse(Tol_Sym.pair_TP4.res$p.value<0.05, "*",
                                            ifelse(Tol_Sym.pair_TP4.res$p.value<0.1, "-", NA))))
Tol_Sym.pair_TP4.res

##Summary statistics by Site and Origin
TP4_TolSym.sum<-summarySE(Therm_TP4, measurevar="Sym.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP4_TolSym.plot<-ggplot(TP4_TolSym.sum, aes(x=Site, y=Sym.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Symbiont Retention", colour="Origin")+
  ylim(0, 1); TP4_TolSym.plot

Chlorophyll

#Check normality
hist(Therm_TP4$Chl.prop)

shapiro.test(Therm_TP4$Chl.prop)

    Shapiro-Wilk normality test

data:  Therm_TP4$Chl.prop
W = 0.94204, p-value = 0.02131

#Not normal

hist(log(Therm_TP4$Chl.prop+1))

shapiro.test(log(Therm_TP4$Chl.prop+1))

    Shapiro-Wilk normality test

data:  log(Therm_TP4$Chl.prop + 1)
W = 0.9512, p-value = 0.04824

##Nearly normal

##Model with log +1 transformation
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_Chl.lme_TP4<-lmer(log(Chl.prop+1)~Origin*Site+(1|Genotype), data=Therm_TP4)

boundary (singular) fit: see help('isSingular')

##Check residuals
Tol_Chl.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_Chl.lme_TP4, plot = F)
plot(Tol_Chl.lme_res_TP4)

##Model results
summary(Tol_Chl.lme_TP4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Chl.prop + 1) ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP4

REML criterion at convergence: -139.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6397 -0.7647 -0.1492  0.5006  2.4251 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.000000 0.00000 
 Residual             0.001803 0.04246 
Number of obs: 47, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              0.09448    0.01226 43.00000   7.707 1.24e-09 ***
OriginTransplant        -0.01266    0.01734 43.00000  -0.730    0.469    
SiteSS                  -0.02659    0.01773 43.00000  -1.500    0.141    
OriginTransplant:SiteSS  0.02939    0.02479 43.00000   1.185    0.242    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.707              
SiteSS      -0.692  0.489       
OrgnTrn:SSS  0.494 -0.699 -0.715
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

eta_squared(Tol_Chl.lme_TP4)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       6.26e-04 | [0.00, 1.00]
Site        |           0.02 | [0.00, 1.00]
Origin:Site |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_TP4.res<-data.frame(summary(Tol_Chl.lme_TP4)$coefficients[-1,])
Tol_Chl_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_Chl_TP4.res$EtaSq<-c(eta_squared(Tol_Chl.lme_TP4)$Eta2)
Tol_Chl_TP4.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl_TP4.res))

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_Chl.pair_TP4<-emmeans(Tol_Chl.lme_TP4, pairwise~Origin | Site)
Tol_Chl.pair_TP4

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.0945 0.0123 21.6   0.0690   0.1199
 Transplant 0.0818 0.0123 21.6   0.0564   0.1073

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native     0.0679 0.0129 23.5   0.0413   0.0944
 Transplant 0.0846 0.0123 21.6   0.0592   0.1101

Degrees-of-freedom method: kenward-roger 
Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE   df t.ratio p.value
 Native - Transplant   0.0127 0.0173 41.0   0.730  0.4693

Site = SS:
 contrast            estimate     SE   df t.ratio p.value
 Native - Transplant  -0.0167 0.0178 41.2  -0.942  0.3517

Note: contrasts are still on the log(mu + 1) scale 
Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_Chl.pair.es_TP4<-data.frame(eff_size(Tol_Chl.pair_TP4, sigma=sigma(Tol_Chl.lme_TP4), edf=df.residual(Tol_Chl.lme_TP4)))

Since 'object' is a list, we are using the contrasts already present.

Tol_Chl.pair.es_TP4

Tol_Chl.pair.es_TP4<-Tol_Chl.pair.es_TP4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_Chl.pair_TP4.res<-merge(data.frame(Tol_Chl.pair_TP4$contrasts), Tol_Chl.pair.es_TP4[,-c(1)])
Tol_Chl.pair_TP4.res$Response<-rep("Chlorophyll Retention", nrow(Tol_Chl.pair_TP4.res))

##Add Significance levels
Tol_Chl.pair_TP4.res$Sig<-ifelse(Tol_Chl.pair_TP4.res$p.value<0.001, "***",
                              ifelse(Tol_Chl.pair_TP4.res$p.value<0.01, "**",
                                     ifelse(Tol_Chl.pair_TP4.res$p.value<0.05, "*",
                                            ifelse(Tol_Chl.pair_TP4.res$p.value<0.1, "-", NA))))
Tol_Chl.pair_TP4.res

##Summary statistics by Site and Origin
TP4_TolChl.sum<-summarySE(Therm_TP4, measurevar="Chl.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP4_TolChl.plot<-ggplot(TP4_TolChl.sum, aes(x=Site, y=Chl.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Chlorophyll Retention", colour="Origin")+
  ylim(0, 0.35); TP4_TolChl.plot

FvFm

#Check normality
hist(Therm_TP4$PAM.prop)

shapiro.test(Therm_TP4$PAM.prop)

    Shapiro-Wilk normality test

data:  Therm_TP4$PAM.prop
W = 0.86379, p-value = 5.159e-05

#Not normal

hist(Therm_TP4$PAM.prop^2)

shapiro.test(Therm_TP4$PAM.prop^2)

    Shapiro-Wilk normality test

data:  Therm_TP4$PAM.prop^2
W = 0.90996, p-value = 0.00134

##Still not normal

##Model
#Function of Site and Origin, with Genotype as a Random effect
#Interaction between Site and Origin
Tol_PAM.lme_TP4<-lmer(PAM.prop^2~Origin*Site+(1|Genotype), data=Therm_TP4)

##Check residuals
Tol_PAM.lme_res_TP4 <- simulateResiduals(fittedModel = Tol_PAM.lme_TP4, plot = F)
plot(Tol_PAM.lme_res_TP4)

##Model results
summary(Tol_PAM.lme_TP4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: PAM.prop^2 ~ Origin * Site + (1 | Genotype)
   Data: Therm_TP4

REML criterion at convergence: -31

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.81938 -0.46352 -0.06318  0.64045  2.58461 

Random effects:
 Groups   Name        Variance Std.Dev.
 Genotype (Intercept) 0.01982  0.1408  
 Residual             0.02031  0.1425  
Number of obs: 48, groups:  Genotype, 3

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)              0.58067    0.09109  2.78330   6.374  0.00976 **
OriginTransplant         0.03985    0.05818 42.00000   0.685  0.49718   
SiteSS                   0.03545    0.05818 42.00000   0.609  0.54566   
OriginTransplant:SiteSS -0.07357    0.08228 42.00000  -0.894  0.37633   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) OrgnTr SiteSS
OrgnTrnspln -0.319              
SiteSS      -0.319  0.500       
OrgnTrn:SSS  0.226 -0.707 -0.707

eta_squared(Tol_PAM.lme_TP4)

# Effect Size for ANOVA (Type III)

Parameter   | Eta2 (partial) |       95% CI
-------------------------------------------
Origin      |       1.32e-04 | [0.00, 1.00]
Site        |       2.53e-05 | [0.00, 1.00]
Origin:Site |           0.02 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_PAM_TP4.res<-data.frame(summary(Tol_PAM.lme_TP4)$coefficients[-1,])
Tol_PAM_TP4.res$Predictor<-c("Origin", "Site", "Origin x Site")
Tol_PAM_TP4.res$EtaSq<-c(eta_squared(Tol_PAM.lme_TP4)$Eta2)
Tol_PAM_TP4.res$Response<-rep("FvFm Retention", nrow(Tol_PAM_TP4.res))

##Pairwise Comparisons of Interest
#Native vs Transplant within a Site
Tol_PAM.pair_TP4<-emmeans(Tol_PAM.lme_TP4, pairwise~Origin | Site)
Tol_PAM.pair_TP4

$emmeans
Site = KL:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.581 0.0911 2.78    0.278    0.884
 Transplant  0.621 0.0911 2.78    0.317    0.924

Site = SS:
 Origin     emmean     SE   df lower.CL upper.CL
 Native      0.616 0.0911 2.78    0.313    0.919
 Transplant  0.582 0.0911 2.78    0.279    0.885

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant  -0.0398 0.0582 42  -0.685  0.4972

Site = SS:
 contrast            estimate     SE df t.ratio p.value
 Native - Transplant   0.0337 0.0582 42   0.580  0.5652

Degrees-of-freedom method: kenward-roger 

#Calculate standardized effect sizes for pairwise comparisons
Tol_PAM.pair.es_TP4<-data.frame(eff_size(Tol_PAM.pair_TP4, sigma=sigma(Tol_PAM.lme_TP4), edf=df.residual(Tol_PAM.lme_TP4)))

Since 'object' is a list, we are using the contrasts already present.

Tol_PAM.pair.es_TP4

Tol_PAM.pair.es_TP4<-Tol_PAM.pair.es_TP4 %>% dplyr::rename(contrast.se = contrast, SE.es = SE, df.es = df)

##Save Results
Tol_PAM.pair_TP4.res<-merge(data.frame(Tol_PAM.pair_TP4$contrasts), Tol_PAM.pair.es_TP4[,-c(1)])
Tol_PAM.pair_TP4.res$Response<-rep("FvFm Retention", nrow(Tol_PAM.pair_TP4.res))

##Add Significance levels
Tol_PAM.pair_TP4.res$Sig<-ifelse(Tol_PAM.pair_TP4.res$p.value<0.001, "***",
                              ifelse(Tol_PAM.pair_TP4.res$p.value<0.01, "**",
                                     ifelse(Tol_PAM.pair_TP4.res$p.value<0.05, "*",
                                            ifelse(Tol_PAM.pair_TP4.res$p.value<0.1, "-", NA))))
Tol_PAM.pair_TP4.res

##Summary statistics by Site and Origin
TP4_TolPAM.sum<-summarySE(Therm_TP4, measurevar="PAM.prop", groupvars=c("Site", "Origin", "Site.Orig"), na.rm=TRUE)

##Plot Average Retention across Treatments
TP4_TolPAM.plot<-ggplot(TP4_TolPAM.sum, aes(x=Site, y=PAM.prop, colour=Site.Orig)) + 
  scale_colour_manual(values=Orig.colors.o)+
  geom_errorbar(aes(ymin=PAM.prop-se, ymax=PAM.prop+se), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(size=point.sz, position=position_dodge(width=0.5))+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x="Site and Origin", y="Photochemical Efficiency Retention", colour="Origin")+
  ylim(0.6, 1); TP4_TolPAM.plot

Save Results

##Combine Results
Tol_TP4.res<-rbind(Tol_PAM_TP4.res, Tol_Sym_TP4.res, Tol_Chl_TP4.res)
Tol.pair.res_TP4<-rbind(Tol_PAM.pair_TP4.res, Tol_Sym.pair_TP4.res, Tol_Chl.pair_TP4.res)

##Add Timepoint
Tol_TP4.res$TimeP<-rep("TP4", nrow(Tol_TP4.res))
Tol.pair.res_TP4$TimeP<-rep("TP4", nrow(Tol.pair.res_TP4))

Percent Difference in Fv/Fm

Focus on Fv/Fm at KL TP 1 and 2 where difference were significant

##Summary statistics by Genotype, Origin, Site, and Timepoint
Tol.PAM.sum<-summarySE(Therm_H, measurevar="PAM.prop", groupvars=c("Genotype", "Origin", "Site", "TimeP"), na.rm=TRUE)

##Remove Initial Timepoint
Tol.PAM.sum<-Tol.PAM.sum[-c(which(Tol.PAM.sum$TimeP=="IN")),]

##Calculate Percent Difference between Native and Transplant by Genotype
Tol.PAM.dif<-Tol.PAM.sum[which(Tol.PAM.sum$Origin=="Native"), c(1,3,4,6)]
Tol.PAM.dif<-Tol.PAM.dif %>% dplyr::rename(PAM_N = PAM.prop)
Tol.PAM.dif$PAM_T<-Tol.PAM.sum$PAM.prop[which(Tol.PAM.sum$Origin=="Transplant")]

##KL
Tol.PAM.dif_KL<-subset(Tol.PAM.dif, Site=="KL")
Tol.PAM.dif_KL$Perc.Inc<-((Tol.PAM.dif_KL$PAM_T-Tol.PAM.dif_KL$PAM_N)/Tol.PAM.dif_KL$PAM_N)*100
Tol.PAM.dif_KL

##Mean and Standard Error by Timepoint
Tol.PAM.dif_KL.sum<-summarySE(Tol.PAM.dif_KL, measurevar="Perc.Inc", groupvars=c("TimeP"), na.rm=TRUE)
Tol.PAM.dif_KL.sum

Performance Effect Sizes

Current vs Original Site

##Combine Results
Growth.TP.res<-rbind(TLE_TP1.2.res, TLE_TP3.4.res)
Growth.TP.res<-Growth.TP.res %>% dplyr::rename(TimeP = Timepoint)
Tol_TP.res<-rbind(Tol_TP1.res, Tol_TP2.res, Tol_TP3.res, Tol_TP4.res)
Perf.res<-rbind(Growth.TP.res, Tol_TP.res)
Perf.res<-Perf.res %>% dplyr::rename(p = Pr...t..)

##Responses for Figure
Perf.res.plot<-Perf.res[-c(which(Perf.res$Response=="Chlorophyll Retention")),]

##Add Significance levels
Perf.res.plot$Sig<-ifelse(Perf.res.plot$p<0.001, "***", ifelse(Perf.res.plot$p<0.01, "**", 
            ifelse(Perf.res.plot$p<0.05, "*", ifelse(Perf.res.plot$p<0.1, "-", NA))))

##Set Order
Perf.res.plot$Response<-factor(Perf.res.plot$Response, levels=c("Growth", "FvFm Retention", "Symbiont Retention"), ordered=TRUE)

##Match Growth Timepoints
Perf.res.plot$TimeP[which(Perf.res.plot$TimeP=="TP1_2")]<-"TP2"
Perf.res.plot$TimeP[which(Perf.res.plot$TimeP=="TP3_4")]<-"TP4"

Plot Effect Size of Origin and Site

Perf.Site.Orig.ES.plot<-ggplot(Perf.res.plot, aes(fill=Predictor, y=EtaSq, x=TimeP)) + 
    geom_bar(position="stack", stat="identity")+
  facet_wrap(vars(Response))+
  theme_bw()+
  scale_fill_manual(values=c("indianred2", "mediumpurple3", "steelblue2"))+
  theme(strip.text=element_text(size=axis.title.sz),
        strip.background = element_rect(fill="white"),
        axis.title.x = element_text(size = axis.title.sz), 
        axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"),
        legend.text=element_text(size=leg.txt.sz), 
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"), 
        legend.position=c(.08, .8), 
        panel.grid.major.x=element_blank())+
   labs(x="Time Point", y=expression(paste("Predictor Effect Size (p", eta^2, ")")), colour="Predictor")+
   geom_text(aes(label=Sig), position=position_stack(vjust=0.5), size=sig.sz); Perf.Site.Orig.ES.plot

Native vs Transplant

##Combine Pairwise Results
Growth.pair.res<-rbind(TLE.pair_TP1.2.res, TLE.pair_TP3.4.res)
Tol.pair.res<-rbind(Tol.pair.res_TP1, Tol.pair.res_TP2, Tol.pair.res_TP3, Tol.pair.res_TP4)
Perf.pair.res<-rbind(Growth.pair.res[,c(1:13, 15, 14)], Tol.pair.res)

##Responses for Figure
Perf.pair.res.plot<-Perf.pair.res[-c(which(Perf.pair.res$Response=="Chlorophyll Retention")),]

##Add Significance
Perf.pair.res.plot$Significant<-ifelse(Perf.pair.res.plot$p.value<0.05, "sig", "non")

##Set Order
Perf.pair.res.plot$Response<-factor(Perf.pair.res.plot$Response, levels=c("Growth", "FvFm Retention", "Symbiont Retention"), ordered=TRUE)

##Match Growth Timepoints
Perf.pair.res.plot$TimeP[which(Perf.pair.res.plot$TimeP=="TP1_2")]<-"TP2"
Perf.pair.res.plot$TimeP[which(Perf.pair.res.plot$TimeP=="TP3_4")]<-"TP4"

##Set Factors
Perf.pair.res.plot$Site<-factor(Perf.pair.res.plot$Site, levels=c("KL", "SS"))

Plot Effect Size of Origin

##Plot Effect Size across Timepoints
Perf.Orig.ES.plot<-ggplot(Perf.pair.res.plot, aes(x=TimeP, y=effect.size)) + 
  geom_rect(fill="grey90", alpha=0.4, xmin=-Inf, xmax=Inf, ymin=-.4, ymax=0.4)+
  geom_rect(fill="grey80", alpha=0.4, xmin=-Inf, xmax=Inf, ymin=-.25, ymax=0.25)+
  geom_rect(fill="grey70", alpha=0.4, xmin=-Inf, xmax=Inf, ymin=-.1, ymax=0.1)+
  geom_hline(yintercept=0, linetype="dashed", color="grey50", linewidth=0.5)+
  geom_errorbar(aes(ymin=effect.size-SE.es, ymax=effect.size+SE.es, colour=Site), width=cap.sz, linewidth=bar.sz, position=position_dodge(width=0.5)) +
  geom_point(aes(shape=Significant, colour=Site), size=point.sz, position=position_dodge(width=0.5))+
  facet_wrap(vars(Response))+
  scale_colour_manual(values=Site.colors.o)+
  scale_shape_manual(values=c(1, 16))+
  theme_bw()+
  theme(axis.title.x = element_text(size = axis.title.sz), 
        axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), 
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"), 
        legend.position=c(.1, .1),
        legend.direction = "horizontal",
        strip.text=element_text(size=axis.title.sz),
        strip.background = element_rect(fill="white"), 
        panel.grid.major.y=element_blank())+
  labs(x="Time Point", y="Origin Effect Size (Cohen's d)", colour="Site")+
  guides(shape="none"); Perf.Orig.ES.plot

Need to fix dodge for points and error bars

Trade Offs

Initial T1 and T2

Average Performance

Calculating the average of Growth (TP1 to TP2) and Thermal Tolerance (Retention of Symbionts, Chlorophyll, and Photochemical Efficiency at TP1) for each Set (Site, Genotype, Origin).

Growth

##Average
names(Growth_TP1.2)

 [1] "Site"       "ID"         "Orig"       "Origin"     "Genotype"   "TL.TP1_cm"  "TL.TP2_cm" 
 [8] "TL.TP3_cm"  "TL.TP4_cm"  "Set"        "Site.Orig"  "TP1"        "TP2"        "TP1.2_days"
[15] "Ext_cm"     "TLE_cm.day"

TP1.2_Growth_Set<-aggregate(Growth_TP1.2$TLE_cm.day, list(Growth_TP1.2$Site, Growth_TP1.2$Genotype, Growth_TP1.2$Orig), mean)
names(TP1.2_Growth_Set)<-c("Site", "Genotype", "Orig", "TLE_cm.day")

##Add Standard Error
TP1.2_Growth_Set$TLE_se<-aggregate(Growth_TP1.2$TLE_cm.day, list(Growth_TP1.2$Site, Growth_TP1.2$Genotype, Growth_TP1.2$Orig), std.error)[,4]

Tolerance

##Average
names(Therm_TP1)

 [1] "TimeP"         "Site"          "Genotype"      "Origin"        "ID"           
 [6] "RandN"         "Treat"         "Treatment"     "Orig"          "Set"          
[11] "Site.Orig"     "SA_cm2"        "Chl_ug.cm2"    "Sym10.6_cm2"   "Fv_Fm"        
[16] "Chl_ug.cm2_C"  "Sym10.6_cm2_C" "Fv_Fm_C"       "Chl.prop"      "Sym.prop"     
[21] "PAM.prop"     

T1_Therm_Set<-aggregate(Therm_TP1[,c(19:21)], list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), mean)
names(T1_Therm_Set)<-c("Site", "Genotype", "Orig", "Chl.prop", "Sym.prop", "PAM.prop")

##Add Standard Error
T1_Therm_Set$Chl_se<-aggregate(Therm_TP1$Chl.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]

T1_Therm_Set$Sym_se<-aggregate(Therm_TP1$Sym.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]

T1_Therm_Set$PAM_se<-aggregate(Therm_TP1$PAM.prop, list(Therm_TP1$Site, Therm_TP1$Genotype, Therm_TP1$Orig), std.error)[,4]

Merge Growth and Tolerance

TP1.2_Perf_Set<-merge(TP1.2_Growth_Set, T1_Therm_Set, all=TRUE)
TP1.2_Perf_Set$Set<-paste(TP1.2_Perf_Set$Site, TP1.2_Perf_Set$Genotype, TP1.2_Perf_Set$Orig, sep=".")

##Set factor variables
TP1.2_Perf_Set$Site<-factor(TP1.2_Perf_Set$Site, levels=c("KL", "SS"))
TP1.2_Perf_Set$Genotype<-factor(TP1.2_Perf_Set$Genotype, levels=c("AC8", "AC10", "AC12"))
TP1.2_Perf_Set$Orig<-factor(TP1.2_Perf_Set$Orig, levels=c("N", "T"))

##Add a Sample Set Variable
TP1.2_Perf_Set$Set<-paste(TP1.2_Perf_Set$Site, TP1.2_Perf_Set$Genotype, TP1.2_Perf_Set$Orig, sep=".")
TP1.2_Perf_Set$Set<-factor(TP1.2_Perf_Set$Set)

##Add Site.Orig variable
TP1.2_Perf_Set$Site.Orig<-paste(TP1.2_Perf_Set$Site, TP1.2_Perf_Set$Orig, sep=".")
TP1.2_Perf_Set$Site.Orig<-factor(TP1.2_Perf_Set$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))

Correlation between Growth and Tolerance

Symbiont Retention

cor.test(TP1.2_Perf_Set$TLE_cm.day, TP1.2_Perf_Set$Sym.prop, method="spearman")

    Spearman's rank correlation rho

data:  TP1.2_Perf_Set$TLE_cm.day and TP1.2_Perf_Set$Sym.prop
S = 462, p-value = 0.03733
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.6153846 

TP1.2_TLE_Sym.plot<-ggplot(TP1.2_Perf_Set, aes(x=TLE_cm.day, y=Sym.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Symbiont Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = -.15, label=expression(bolditalic(paste(r[S], " = -0.615, p = 0.037"))), size=sig.sz, hjust = 0); TP1.2_TLE_Sym.plot

Chlorophyll Retention

cor.test(TP1.2_Perf_Set$TLE_cm.day, TP1.2_Perf_Set$Chl.prop, method="spearman")

    Spearman's rank correlation rho

data:  TP1.2_Perf_Set$TLE_cm.day and TP1.2_Perf_Set$Chl.prop
S = 454, p-value = 0.04884
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.5874126 

TP1.2_TLE_Chl.plot<-ggplot(TP1.2_Perf_Set, aes(x=TLE_cm.day, y=Chl.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Chlorophyll Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = -.05, label=expression(bolditalic(paste(r[S], " = -0.587, p = 0.049"))), size=sig.sz, hjust = 0); TP1.2_TLE_Chl.plot

FvFm Retention

cor.test(TP1.2_Perf_Set$TLE_cm.day, TP1.2_Perf_Set$PAM.prop, method="spearman")

    Spearman's rank correlation rho

data:  TP1.2_Perf_Set$TLE_cm.day and TP1.2_Perf_Set$PAM.prop
S = 480, p-value = 0.01883
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.6783217 

TP1.2_TLE_PAM.plot<-ggplot(TP1.2_Perf_Set, aes(x=TLE_cm.day, y=PAM.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Photo. Effic. Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.7, label=expression(bolditalic(paste(r[S], " = -0.678, p = 0.019"))), size=sig.sz, hjust = 0); TP1.2_TLE_PAM.plot

Final T3 and T4

Average Performance

Calculating the average of Growth (TP3 to TP4) and Thermal Tolerance (Retention of Symbionts, Chlorophyll, and Photochemical Efficiency at TP3) for each Set (Site, Genotype, Origin).

Growth

##Average
names(Growth_TP3.4)

 [1] "Site"       "ID"         "Orig"       "Origin"     "Genotype"   "TL.TP1_cm"  "TL.TP2_cm" 
 [8] "TL.TP3_cm"  "TL.TP4_cm"  "Set"        "Site.Orig"  "TP3"        "TP4"        "TP3.4_days"
[15] "Ext_cm"     "TLE_cm.day"

TP3.4_Growth_Set<-aggregate(Growth_TP3.4$TLE_cm.day, list(Growth_TP3.4$Site, Growth_TP3.4$Genotype, Growth_TP3.4$Orig), mean)
names(TP3.4_Growth_Set)<-c("Site", "Genotype", "Orig", "TLE_cm.day")

##Add Standard Error
TP3.4_Growth_Set$TLE_se<-aggregate(Growth_TP3.4$TLE_cm.day, list(Growth_TP3.4$Site, Growth_TP3.4$Genotype, Growth_TP3.4$Orig), std.error)[,4]

Tolerance

##Average
names(Therm_TP4)

 [1] "TimeP"         "Site"          "Genotype"      "Origin"        "ID"           
 [6] "RandN"         "Treat"         "Treatment"     "Orig"          "Set"          
[11] "Site.Orig"     "SA_cm2"        "Chl_ug.cm2"    "Sym10.6_cm2"   "Fv_Fm"        
[16] "Chl_ug.cm2_C"  "Sym10.6_cm2_C" "Fv_Fm_C"       "Chl.prop"      "Sym.prop"     
[21] "PAM.prop"     

T4_Therm_Set<-aggregate(Therm_TP4[,c(19:21)], list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), mean)
names(T4_Therm_Set)<-c("Site", "Genotype", "Orig", "Chl.prop", "Sym.prop", "PAM.prop")

##Add Standard Error
T4_Therm_Set$Chl_se<-aggregate(Therm_TP4$Chl.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]

T4_Therm_Set$Sym_se<-aggregate(Therm_TP4$Sym.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]

T4_Therm_Set$PAM_se<-aggregate(Therm_TP4$PAM.prop, list(Therm_TP4$Site, Therm_TP4$Genotype, Therm_TP4$Orig), std.error)[,4]

Merge Growth and Tolerance

TP3.4_Perf_Set<-merge(TP3.4_Growth_Set, T4_Therm_Set, all=TRUE)
TP3.4_Perf_Set$Set<-paste(TP3.4_Perf_Set$Site, TP3.4_Perf_Set$Genotype, TP3.4_Perf_Set$Orig, sep=".")

##Set factor variables
TP3.4_Perf_Set$Site<-factor(TP3.4_Perf_Set$Site, levels=c("KL", "SS"))
TP3.4_Perf_Set$Genotype<-factor(TP3.4_Perf_Set$Genotype, levels=c("AC8", "AC10", "AC12"))
TP3.4_Perf_Set$Orig<-factor(TP3.4_Perf_Set$Orig, levels=c("N", "T"))

##Add a Sample Set Variable
TP3.4_Perf_Set$Set<-paste(TP3.4_Perf_Set$Site, TP3.4_Perf_Set$Genotype, TP3.4_Perf_Set$Orig, sep=".")
TP3.4_Perf_Set$Set<-factor(TP3.4_Perf_Set$Set)

##Add Site.Orig variable
TP3.4_Perf_Set$Site.Orig<-paste(TP3.4_Perf_Set$Site, TP3.4_Perf_Set$Orig, sep=".")
TP3.4_Perf_Set$Site.Orig<-factor(TP3.4_Perf_Set$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))

Correlation between Growth and Tolerance

Symbiont Retention

cor.test(TP3.4_Perf_Set$TLE_cm.day, TP3.4_Perf_Set$Sym.prop, method="spearman")

    Spearman's rank correlation rho

data:  TP3.4_Perf_Set$TLE_cm.day and TP3.4_Perf_Set$Sym.prop
S = 444, p-value = 0.06663
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.5524476 

TP3.4_TLE_Sym.plot<-ggplot(TP3.4_Perf_Set, aes(x=TLE_cm.day, y=Sym.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Symbiont Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.1, label=expression(italic(paste(r[S], " = -0.552, p = 0.067"))), size=sig.sz, hjust = 0); TP3.4_TLE_Sym.plot

Chlorophyll Retention

cor.test(TP3.4_Perf_Set$TLE_cm.day, TP3.4_Perf_Set$Chl.prop, method="spearman")

    Spearman's rank correlation rho

data:  TP3.4_Perf_Set$TLE_cm.day and TP3.4_Perf_Set$Chl.prop
S = 148, p-value = 0.327
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.3272727 

TP3.4_TLE_Chl.plot<-ggplot(TP3.4_Perf_Set, aes(x=TLE_cm.day, y=Chl.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Chlorophyll Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.065, label=expression(italic(paste(r[S], " = 0.327, p = 0.327"))), size=sig.sz, hjust = 0); TP3.4_TLE_Chl.plot

FvFm Retention

cor.test(TP3.4_Perf_Set$TLE_cm.day, TP3.4_Perf_Set$PAM.prop, method="spearman")

    Spearman's rank correlation rho

data:  TP3.4_Perf_Set$TLE_cm.day and TP3.4_Perf_Set$PAM.prop
S = 500, p-value = 0.007353
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.7482517 

TP3.4_TLE_PAM.plot<-ggplot(TP3.4_Perf_Set, aes(x=TLE_cm.day, y=PAM.prop)) + 
   geom_smooth(method="lm", color="#AA185AFF", fill="#F8D7BFFF", se=TRUE)+
   geom_point(aes(colour=Site.Orig, shape=Genotype),size=point.sz)+
    scale_colour_manual(values=Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
   theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"), axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz), legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"))+
  labs(x=expression(paste('Growth Rate (cm day'^-1*")")), 
       y="Photo. Effic. Retention", 
       colour="Site Origin")+
   annotate("text", x = 0.07, y = 0.5, label=expression(bolditalic(paste(r[S], " = -0.748, p = 0.007"))), size=sig.sz, hjust = 0); TP3.4_TLE_PAM.plot

Figures

Figure 4 Performance Effect Sizes

Adjust for Panel

Perf.Site.Orig.ES.plot<-Perf.Site.Orig.ES.plot+ labs(x="")

Figure 4

##Create Panel
Perf.ES_fig<-plot_grid(Perf.Site.Orig.ES.plot, Perf.Orig.ES.plot, 
                    rel_widths=1, rel_heights = 1,
                    nrow=2, ncol=1, byrow=T, labels = c("A", "B"), 
                    label_size=panel.lab.sz)

Warning: Removed 18 rows containing missing values (`geom_text()`).

##Save Figure
ggsave(filename="Figures/03_Performance/Fig4_PerformanceEffectSizes.png", plot=Perf.ES_fig, dpi=300, width=12, height=9, units="in")

Figure 5 Growth vs Tolerance

Adjust Figures for Panel

##Timepoint 1 to 2
TP1.2_TLE_PAM.plot<-TP1.2_TLE_PAM.plot + ggtitle("Photo. Effic.")+ labs(x=NULL, y = "Tolerance (Retention)")+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5), legend.position="none")

TP1.2_TLE_Sym.plot<- TP1.2_TLE_Sym.plot + ggtitle("Symbionts")+ labs(x=NULL, y = NULL)+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, hjust=0.5)) + ylim(-.15, .8)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

##Timepoint 3 to 4
TP3.4_TLE_PAM.plot<-TP3.4_TLE_PAM.plot + labs(y = "Tolerance (Retention)") +
  theme(legend.position="none")

TP3.4_TLE_Sym.plot<-TP3.4_TLE_Sym.plot + labs(y = NULL)

Figure 5

##Create Panel
TradeOff_fig<-plot_grid(TP1.2_TLE_PAM.plot, TP1.2_TLE_Sym.plot,
                      TP3.4_TLE_PAM.plot, TP3.4_TLE_Sym.plot,
                    rel_widths=c(.75, 1), rel_heights = c(1, .98),
                    nrow=2, ncol=2, byrow=T, 
                    labels = c("A", "B", "C", "D"), 
                    label_size=panel.lab.sz)

`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).Warning: Removed 1 rows containing missing values (`geom_point()`).Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'`geom_smooth()` using formula = 'y ~ x'Warning: is.na() applied to non-(list or vector) of type 'expression'

##Save Figure
ggsave(filename="Figures/03_Performance/Fig5_TradeOff.png", plot=TradeOff_fig, dpi=300, width=8, height=8, units="in")

Figure S5 Growth

Adjust Figures for Panel

##TP 1 to TP 2
TP1.2_Growth.plot<-TP1.2_Growth.plot+ ggtitle("TP1 to TP2")+
  theme(plot.title = element_text(colour="black", 
                                  size=panel.lab.sz, hjust=0.5), 
        legend.position=c(0.25, 0.85), legend.direction="horizontal")+
  labs(colour=NULL)+
  guides(color=guide_legend(nrow=2, byrow=FALSE)) 


##TP 3 to TP 4
TP3.4_Growth.plot<-TP3.4_Growth.plot+ ggtitle("TP3 to TP4")+
  theme(plot.title = element_text(colour="black", 
                                  size=panel.lab.sz, hjust=0.5), 
        legend.position="none")

Figure S5

##Create Panel
Growth_fig<-plot_grid(TP1.2_Growth.plot, TP3.4_Growth.plot, 
                    rel_widths=1, rel_heights = 1,
                    nrow=1, ncol=2, byrow=T, labels = c("A", "B"), 
                    label_size=panel.lab.sz)

##Save Figure
ggsave(filename="Figures/03_Performance/FigS5_Growth.png", plot=Growth_fig, dpi=300, width=10, height=4.5, units="in")

Figure S6 Initial Heat Assay

Adjust Figures for Panel

##FvFm
IN_PAM.plot<-IN_PAM.plot + labs(x="", y="Fv/Fm") + theme(legend.position="none") + ylim(0.5, 0.655)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

IN_TolPAM.plot<-IN_TolPAM.plot + labs(x="", y="Retention") + 
  theme(legend.position="none") + ylim(0.75, 1)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

##Symbionts
IN_Sym.plot<-IN_Sym.plot + labs(x="", y="Symbionts") + theme(legend.position="none") + ylim(0,1)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

IN_TolSym.plot<-IN_TolSym.plot + labs(x="", y="Retention") + 
  theme(legend.position="none") + ylim(.25, .75)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

##Chlorophyll
IN_Chl.plot<-IN_Chl.plot + labs(x="Site and Treatment", y="Chlorophyll") + 
  theme(legend.position="none") + ylim(0,1.25)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

IN_TolChl.plot<-IN_TolChl.plot + labs(x="Site", y="Retention") + 
  theme(legend.position="none") + ylim(0, 0.25)

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Figure S6

##Create Panel
IN_HeatAssay_fig<-plot_grid(IN_PAM.plot, IN_TolPAM.plot, 
                      IN_Sym.plot, IN_TolSym.plot,
                      IN_Chl.plot, IN_TolChl.plot,
                    rel_widths=c(1, 0.75), 
                    rel_heights = c(1, 1, 1),
                    nrow=3, ncol=2, byrow=T, 
                    labels = c("A", "B", "C", "D", "E", "F"), 
                    label_size=panel.lab.sz)

##Save Figure
ggsave(filename="Figures/03_Performance/FigS6_InitialHeatAssay.png", plot=IN_HeatAssay_fig, dpi=300, width=8, height=11, units="in")

Figure S7 Thermal Tolerance

Adjust Figures for Panel

##FvFm Retention
TP1_TolPAM.plot<-TP1_TolPAM.plot + ggtitle("TP1")+ 
  labs(x=NULL, y="Photo. Effic. Retention", colour=NULL)+
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position=c(0.5, 0.15), legend.direction="horizontal")+
  guides(color=guide_legend(nrow=2, byrow=FALSE))

TP2_TolPAM.plot<-TP2_TolPAM.plot+ ggtitle("TP2")+ labs(x=NULL, y=NULL) +
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position="none")

TP3_TolPAM.plot<-TP3_TolPAM.plot+ ggtitle("TP3")+ labs(x=NULL, y=NULL) +
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position="none")

TP4_TolPAM.plot<-TP4_TolPAM.plot+ ggtitle("TP4")+ labs(x=NULL, y=NULL) +
  theme(plot.title = element_text(colour="black", size=panel.lab.sz, face="bold", hjust=0.5), 
        legend.position="none")


##Symbiont Retention
TP1_TolSym.plot<- TP1_TolSym.plot + labs(x=NULL) + theme(legend.position="none")

TP2_TolSym.plot<- TP2_TolSym.plot + labs(x=NULL, y=NULL) + theme(legend.position="none")

TP3_TolSym.plot<- TP3_TolSym.plot + labs(x=NULL, y=NULL) + theme(legend.position="none")

TP4_TolSym.plot<- TP4_TolSym.plot + labs(x=NULL, y=NULL) + theme(legend.position="none")


##Chlorophyll Retention
TP1_TolChl.plot<-TP1_TolChl.plot + theme(legend.position="none")

TP2_TolChl.plot<-TP2_TolChl.plot + labs(y=NULL) + theme(legend.position="none")

TP3_TolChl.plot<-TP3_TolChl.plot + labs(y=NULL) + theme(legend.position="none")

TP4_TolChl.plot<-TP4_TolChl.plot + labs(y=NULL) + theme(legend.position="none")

Figure S7

##Create Panel
Tolerance_fig<-plot_grid(TP1_TolPAM.plot, TP2_TolPAM.plot, 
                      TP3_TolPAM.plot, TP4_TolPAM.plot, 
                      TP1_TolSym.plot, TP2_TolSym.plot, 
                      TP3_TolSym.plot, TP4_TolSym.plot, 
                      TP1_TolChl.plot, TP2_TolChl.plot, 
                      TP3_TolChl.plot, TP4_TolChl.plot, 
                    rel_widths=1, rel_heights = 1,
                    nrow=3, ncol=4, byrow=T, labels = NULL)

##Save Figure
ggsave(filename="Figures/03_Performance/FigS7_Tolerance.png", plot=Tolerance_fig, dpi=300, width=14, height=12, units="in")

Tables

Table S6 Performance LM

##Combine Results Tables
TableS6_Perf.LM<-Perf.res
  
##Organize
names(TableS6_Perf.LM)

[1] "Estimate"   "Std..Error" "df"         "t.value"    "p"          "Predictor"  "EtaSq"     
[8] "Response"   "TimeP"     

TableS6_Perf.LM<-TableS6_Perf.LM %>% dplyr::rename(SE = Std..Error, DF = df, t = t.value, Timepoint = TimeP)
TableS6_Perf.LM<-TableS6_Perf.LM[,c("Timepoint", "Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EtaSq")]

#Round to 3 digits
TableS6_Perf.LM$Estimate<-round(TableS6_Perf.LM$Estimate, 3)
TableS6_Perf.LM$SE<-round(TableS6_Perf.LM$SE, 3)
TableS6_Perf.LM$t<-round(TableS6_Perf.LM$t, 3)
TableS6_Perf.LM$p<-round(TableS6_Perf.LM$p, 3)
TableS6_Perf.LM$EtaSq<-round(TableS6_Perf.LM$EtaSq, 3)

#Integer
TableS6_Perf.LM$DF<-round(TableS6_Perf.LM$DF, 0)

##Write Out Table
write.csv(TableS6_Perf.LM, "Tables/TableS6_Performance_LM_Results.csv", row.names=FALSE)

Table S7 Performance Pairwise

##Save Results Table
TableS7_Perf.pair<-Perf.pair.res

##Organize
names(TableS7_Perf.pair)

 [1] "Site"        "contrast"    "estimate"    "SE"          "df"          "t.ratio"    
 [7] "p.value"     "effect.size" "SE.es"       "df.es"       "lower.CL"    "upper.CL"   
[13] "Response"    "Sig"         "TimeP"      

TableS7_Perf.pair<-TableS7_Perf.pair %>% dplyr::rename(Timepoint = TimeP, Contrast =contrast, Estimate = estimate, DF = df, t = t.ratio, p = p.value, EffectSize = effect.size, EffectSize.SE = SE.es)
TableS7_Perf.pair<-TableS7_Perf.pair[,c("Timepoint", "Response", "Site", "Contrast", "Estimate", "SE", "DF", "t", "p", "EffectSize", "EffectSize.SE")]

#Round to 3 digits
TableS7_Perf.pair$Estimate<-round(TableS7_Perf.pair$Estimate, 3)
TableS7_Perf.pair$SE<-round(TableS7_Perf.pair$SE, 3)
TableS7_Perf.pair$t<-round(TableS7_Perf.pair$t, 3)
TableS7_Perf.pair$p<-round(TableS7_Perf.pair$p, 3)
TableS7_Perf.pair$EffectSize<-round(TableS7_Perf.pair$EffectSize, 3)
TableS7_Perf.pair$EffectSize.SE<-round(TableS7_Perf.pair$EffectSize.SE, 3)

#Integer
TableS7_Perf.pair$DF<-round(TableS7_Perf.pair$DF, 0)

##Write Out Table
write.csv(TableS7_Perf.pair, "Tables/TableS7_Performance_Pairwise_Results.csv", row.names=FALSE)

Table S8 Initial Heat Assay LM

##Combine Results Tables
TableS8_IN_Tolerance<-data.frame(rbind(PAM_IN.res, Sym_IN.res, Chl_IN.res))

##Organize
names(TableS8_IN_Tolerance)

[1] "estimate"  "SE"        "df"        "t.ratio"   "p.value"   "Predictor" "Response" 
[8] "EtaSq"    

TableS8_IN_Tolerance<-TableS8_IN_Tolerance %>% dplyr::rename(Estimate = estimate, DF = df, t = t.ratio, p = p.value, EffectSize = EtaSq)
TableS8_IN_Tolerance<-TableS8_IN_Tolerance[,c("Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EffectSize")]

#Round to 3 digits
TableS8_IN_Tolerance$Estimate<-round(TableS8_IN_Tolerance$Estimate, 3)
TableS8_IN_Tolerance$SE<-round(TableS8_IN_Tolerance$SE, 3)
TableS8_IN_Tolerance$t<-round(TableS8_IN_Tolerance$t, 3)
TableS8_IN_Tolerance$p<-round(TableS8_IN_Tolerance$p, 3)
TableS8_IN_Tolerance$EffectSize<-round(TableS8_IN_Tolerance$EffectSize, 3)

#Integer
TableS8_IN_Tolerance$DF<-round(TableS8_IN_Tolerance$DF, 0)

##Write Out Table
write.csv(TableS8_IN_Tolerance, "Tables/TableS8_Initial_Tolerance_LM_Results.csv", row.names=FALSE)